text:
# Monoclinic and triclinic phases in higher-order Devonshire theory

## I Introduction

Many of the most important and best studied ferroelectric materials adopt the simple-cubic perovskite structure at high temperature and undergo structural phase transitions to distorted ferroelectric structures at lower temperature. Among the best known simple compounds of this kind are BaTiO3 and PbTiO3. Upon cooling, BaTiO3 undergoes a sequence of ferroelectric transitions, first from the cubic (C) to a tetragonal (T), then to an orthorhombic (O), and finally to a rhombohedral (R) phase. Passing through this sequence, the polarization $`𝐏`$ first vanishes in the C phase, and then becomes oriented in the , , and  directions in the T, O, and R phases, respectively. PbTiO3 undergoes a single ferroelectric transition from the C to T phase.

In a classic 1948 paper, Devonshire  was able to explain the observed phases and phase transition sequence quite naturally in terms of a phenomenological Landau-type expansion of the free energy in terms of the ferroelectric order parameter $`𝐏`$. Making use of cubic symmetry and truncating the expansion to sixth order in $`𝐏`$, Devonshire was able to arrive at a simple model with only a single temperature-dependent second-order coefficient, and only three temperature-independent higher-order coefficients. The polarization $`𝐏`$ is the primary order parameter, and the crystallographic labels (T, R, etc.) refer to the distortions induced by the polarization and the resulting strain. Despite its simplicity, this model could successfully reproduce the phase-transition sequence and the piezoelectric and other properties  of BaTiO3. With a simple modification of the anharmonic coefficients, the qualitative behavior of PbTiO3 could be equally well reproduced.

However, the material that is currently in most widespread use for piezoelectric transducer and related applications is the solid solution PbZr1-xTixO3, commonly known as PZT. The standard understanding of the phase diagram of PZT has been as follows. PZT undergoes a transition from the simple cubic C phase to a ferroelectric phase at a Curie temperature $`T_c`$ that ranges from about 490C at $`x`$=1 to 230C at $`x`$=0. The transition occurs to the T phase for $`x`$ greater than about 0.48, and to the R phase for smaller $`x`$. At $`T<150^{}`$C or $`x<0.1`$, some more complex phases involving antiferroelectric (AFE) or antiferrodistortive (AFD, i.e., involving rotations of oxygen octahedra) displacements also occur. The phase boundary between the simple T and R phases, known as the morphotropic phase boundary (MPB), is an almost vertical line in the $`x`$$`T`$ plane at about $`x`$=0.48. Haun et al. have successfully extended the Devonshire model to the case of PZT by including the AFE and AFD degrees of freedom in the phenomenological free energy, still including only terms up to sixth order overall. This model successfully described the simple R–T transition that was understood to occur below $`T_c`$ at the MPB.

As it happened, a surprise was in store. By working with highly purified and carefully prepared samples of PZT, Noheda and coworkers have recently shown that a sliver of monoclinic (M) phase actually interposes itself between the R and T phases in a very narrow composition range (of order 3-4% in $`x`$). That is, at least below $``$100C, the transition is first from T to M at an $`x_{c2}`$ between about 0.48 and 0.51 (depending on temperature), and then from M to R at $`x_{c1}0.47`$, with decreasing $`x`$. The orientation of $`𝐏`$ is, respectively, along , \[$`uuv`$\] ($`u<v`$), and  in the T, M, and R phases, respectively. The experiments have not yet clarified whether or not a direct T–R transition occurs in the higher temperature range 100C$`<T<T_c`$, or whether the sliver of M phase instead persists up to the Curie temperature $`T_c370^{}`$C. Neither the Devonshire theory, nor the modification of Haun et al., predicted the possible occurrence of the M phase. On the other hand, simulations based on a first-principles effective Hamiltonian approach have very recently provided confirmation of the existence of the M phase in just the observed composition range.

We make a brief aside to establish notation. In the C phase (space group $`Pm3m`$), $`𝐏=\mathrm{𝟎}`$. When $`𝐏`$ is constrained to a symmetry axis lying along , , or , the resulting phase and space-group labels become T ($`P4mm`$), R ($`R3m`$), or O ($`Amm2`$), respectively. Similarly, the M phases arise when $`𝐏`$ is confined to a mirror plane. We can distinguish three cases: MC ($`Pm`$), in which $`𝐏`$ is along \[$`0uv`$\]; and MA and MB (both $`Bm`$, sometimes also denoted $`Cm`$) in which $`𝐏`$ lies along \[$`uuv`$\], with $`u<v`$ and $`u>v`$, respectively. The newly observed phase of PZT is of type MA. If $`𝐏`$ is unconstrained by symmetry, strain coupling leads to a triclinic phase ($`P1`$). This exhausts the possible reduced-symmetry states of a cubic perovskite crystal generated by the emergence a single non-zero ferroelectric order parameter, although other phases can occur if AFE and AFD distortions are also present.

Relatively few examples are known of low-symmetry ferroelectrics in which $`𝐏`$ is only constrained to a symmetry plane, or in which $`𝐏`$ is unconstrained by symmetry. The discovery of the $`M_\mathrm{A}`$ phase in a cubic perovskite is thus of considerable note, even aside from the fact that it had been missed for so long and aside from the potential importance of this phase for understanding the large piezoelectric response in PZT.

The failure of the phenomenological theories of Devonshire  and Haun et al.  to describe the existence of the observed MA phase raises an interesting question: What is the simplest and most natural phenomenological model that does predict such a monoclinic phase? Recent work of Souza Filho et al. confirms that the sixth-order Devonshire expansion does not allow for the occurrence of a monoclinic phase and suggests instead a model in which the vanishing of a shear elastic constant at the critical temperature drives the transition to the MA phase. However, this requires the introduction of an additional instability that is unrelated to the ferroelectric one.  That this is unnecessary is demonstrated by the work of Bellaiche, García and Vanderbilt. In their simulations, the shear moduli are taken as temperature-independent, and yet the transition to the monoclinic phase still occurs, driven by the tilting of the ferroelectric polarization away from the symmetry axis.

Suspecting instead that the problem is simply related to the truncation of the expansions of Devonshire and Haun et al. to sixth order, we consider the addition of terms of eighth and higher order to the Devonshire model. The following questions then arise. At what order in the expansion do monoclinic phases first appear in the phase diagram? And, for that matter, at what order do triclinic phases first appear?

The phase diagram of such a model consists of fields in the parameter space of the model (labeled C, T, O, R, M, etc.) within which the order parameter has the specified symmetry (that is, the absolute minimum of the free energy occurs at a point in the order-parameter space of that symmetry). On a trajectory in the phase space which crosses a phase boundary, that global minimum can change into a local minimum, a saddle point, or a maximum, or simply disappear. Understanding how the set of critical or stationary points of the free energy (its minima, saddle points, and maxima in the order-parameter space) varies with the parameters of the model can thus be an important aid in understanding the phase diagram. So, in addition to searches for the global minima by direct computation, we use here a topological analysis of the Devonshire theory to elucidate the critical-point structure of the free energy and answer the above questions.

We show here that the simplest extension of the Devonshire theory, to eighth order in the ferroelectric order parameter, naturally admits all three kinds (MA, MB, MC) of monoclinic phases. Moreover, the eighth-order model makes specific and testable predictions for the types of phase boundaries that can occur. For example, it predicts that the R–MA and T–MA transitions should be of first and second order, respectively. On the other hand, we find that the model has to be extended all the way to twelfth order in order to describe a triclinic ferroelectric phase.

We therefore suggest that the most natural explanation for the occurrence of the MA phase in PZT is that the free-energy surface is unusually anharmonic in this material, such that the eighth-order terms play an important role. It is not difficult to speculate why this might be the case. First, PZT is a disordered material; averaging over the chemical disorder, which plays the role of a quenched random field, may tend to generate higher orders in the phenomenological energy expansion. Second, as will be discussed in Sec. V, there is considerable evidence that the ferroelectric transitions in PZT have a strong order-disorder character; mapping onto a displacive picture may then also tend to generate higher-order terms. And finally, of course, thermal fluctuations and coupling to strain may play some role. In any case, we show below that once one accepts this simple hypothesis of the importance of eighth-order terms in the free-energy expansion, then the behavior of the monoclinic phase observed in experiments  and simulations  can be understood quite naturally.

This paper is organized as follows. Section II establishes the notation used for the expansion of the free energy, and reviews the symmetry considerations that lead to a simplified form of this expansion. The rules governing the numbers and types of stationary points that may occur in the order-parameter space are reviewed in Sec. III. The behavior of the models obtained by truncating the expansion at higher and higher order (fourth, sixth, eighth, tenth, and twelfth) are then carefully elucidated in Sec. IV. Section V gives a brief discussion of a microscopic model in which the T and R phases can be regarded as arising from fluctuations among neighboring, symmetry-equivalent local MA states. Finally, we conclude with a brief summary and discussion of future prospects in Sec. VI.

## II Formalism

We consider the case of a structural phase transition governed by a single continuous vector order parameter $`𝐮`$, such that the free energy

$$F(𝐮,\sigma _i,T)=E(𝐮,\eta _i,S)\underset{i}{}\sigma _i\eta _iTS$$
(1)
is symmetric with respect to operations of the cubic point group. Here $`\sigma _i`$ and $`\eta _i`$ are the stress and strain tensors in Voigt notation, and $`T`$ and $`S`$ are temperature and entropy, respectively. We have in mind primarily the case in which $`𝐮`$ is the ferroelectric polarization $`𝐏`$ in a member of the cubic perovskite class, such as BaTiO3 or PZT, but the use of the symbol $`𝐮`$ serves as a reminder that the formalism applies to a variety of other cases as well. For a crystal with stress-free boundary conditions at some given temperature, we can expand

$`F(𝐮)`$ $`=C_{000}+C_{200}(u_x^2+u_y^2+u_z^2)+C_{400}(u_x^4+u_y^4+u_z^4)`$ (5)
$`+C_{220}(u_x^2u_y^2+u_x^2u_z^2+u_y^2u_z^2)+C_{600}(u_x^6+u_y^6+u_z^6)`$
$`+C_{420}(u_x^4[u_y^2+u_z^2]+u_y^4[u_x^2+u_z^2]+u_z^4[u_x^2+u_y^2])`$
$`+C_{222}u_x^2u_y^2u_z^2+\mathrm{},`$
where terms up to sixth order in $`u`$ have been written explicitly, and the coefficients have been renormalized to subsume the couplings to strain.

For the ferroelectric phases we simplify matters further by focusing on the energy as a function of the orientation of the vector order parameter and introducing the function

$$G(\widehat{𝐮})=\underset{𝐮\widehat{𝐮}}{\mathrm{min}}F(𝐮).$$
(6)
Thus, $`G(\widehat{𝐮})`$ represents the ground-state energy subject to the constraint that the order parameter has given orientation. In most cases of interest, it is reasonable to expect that $`G(\widehat{𝐮})`$ will be a smooth function of $`\widehat{𝐮}`$. In this case, and suppressing the uninteresting constant term in the expansion, it follows that

$`G(\widehat{𝐮})=`$ $`a_4g_4`$ $`+a_6g_6+a_8g_4^2`$ (7)
$`+`$ $`a_{10}g_4g_6+a_{12}g_4^3+a_{12}^{}g_6^2+\mathrm{},`$ (8)
where

$`g_4(\widehat{𝐮})`$ $`=`$ $`x^2y^2+x^2z^2+y^2z^2`$ (9)
$`g_6(\widehat{𝐮})`$ $`=`$ $`x^2y^2z^2`$ (10)
and $`\widehat{𝐮}=(x,y,z)`$ with $`x^2+y^2+z^2=1`$. All independent, symmetry-allowed terms up to twelfth order are explicitly given in Eq. (8); higher orders will not be needed here.

The particular form of the cubic invariants appearing in Eq. (8) is largely arbitrary. For example, one could use the “kubic harmonics”  instead. In Sec. IV D, we shall make use of the expansion of Eqs. (8-10) above. However, for the presentation of numerical results to be given in Sec. IV C, we find it more convenient instead to use the expansion

$$G(\widehat{𝐮})=c_4f_4+c_6f_6+c_8f_8+\mathrm{},$$
(11)
where

$`f_4(\widehat{𝐮})`$ $`=`$ $`12g_4`$ (12)
$`f_6(\widehat{𝐮})`$ $`=`$ $`4g_436g_6`$ (13)
$`f_8(\widehat{𝐮})`$ $`=`$ $`48g_4^212g_436g_6.`$ (14)
Aside from normalization constants, which are chosen to make the range of each function roughly of order unity, this choice can be uniquely defined by the following requirements. (i) The invariant $`f_n`$ contains no terms of order higher than $`n`$. (ii) All three invariants vanish identically for a “tetragonal” (T) value of the order parameter, e.g., $`\widehat{𝐮}=(100)`$. (iii) $`f_6`$ and $`f_8`$ still vanish for a “rhombohedral” (R) value of the order parameter, e.g., $`(111)/\sqrt{3}`$. (iv) $`f_8`$ alone vanishes for an “orthorhombic” (O) value of the order parameter, e.g., $`(110)/\sqrt{2}`$. The values of these functions evaluated at symmetry directions are summarized in Table I for later reference. Requirement (ii) just reflects an arbitrary choice of zero; requirements (iii) and (iv) simplify some later discussion. For example, phase boundaries at which $`G`$(T), $`G`$(R), and $`G`$(O) become degenerate in pairs are easily located (see Sec. IV C).

A few comments about the transition from $`F(𝐮)`$ to $`G(\widehat{𝐮})`$ in Eq. (6) are in order. First, our “reduced” theory in terms of $`G(\widehat{𝐮})`$ can only describe the transitions among the ferroelectric phases, not the transition to the high-temperature cubic phase ($`𝐮=0`$). Second, the minimization over the magnitude of $`𝐮`$ in Eq. (6) generally introduces higher orders into the reduced theory. For example, truncation of Eq. (5) at sixth order may still lead to terms of eighth and higher order in Eqs. (8) and (11). Thus, a statement about what kinds of critical points can occur simultaneously in an $`n^{\mathrm{th}}`$-order version of the reduced theory does not necessarily carry over to the $`n^{\mathrm{th}}`$-order version of the standard Devonshire theory. Nevertheless, a statement that a certain minimum order of expansion is needed for the existence of a certain kind of stationary point in the reduced theory does carry over to the standard Devonshire theory. For example, we show later that the reduced theory must be carried to eighth order to allow for a monoclinic minimum, i.e., an equilibrium $`𝐮`$ lying in a (1$`\overline{1}`$0) plane. Then the same statement applies to the standard Devonshire theory. For, suppose that $`u_0`$ is the magnitude of $`𝐮`$ at the minimum, and let $`G_0(\widehat{𝐮})=F(u_0\widehat{𝐮})`$; then $`G_0`$ must have a monoclinic minimum. But $`G_0`$ is also of no higher than eighth order, and so admits only the same type of minima as does $`G`$.

## III Topological considerations

The rules governing the numbers and types of critical or stationary points of a scalar function defined on a continuous manifold emerge from a branch of algebraic topology called Morse theory or, equivalently, the calculus of variations in the large. These rules have proved useful in analyzing lattice vibration spectra and the excitation spectra of crystalline materials generally. They have not, however, enjoyed the broad application in condensed matter physics that might have been expected from these early successes.

The rules take the form of a set of inequalities and one equality. The latter is the most powerful, and we confine ourselves to considering only it explicitly. $`G(\widehat{𝐮})`$ is an analytic function defined on the surface of a sphere, a closed two-dimensional manifold of genus zero, for which the equality becomes

$$N_0N_1+N_2=2,$$
(15)
where $`N_0`$, $`N_1`$, and $`N_2`$ are the number of minima, saddle points, and maxima, respectively. This formula assumes only analyticity in $`G(\widehat{𝐮})`$ and that its Hessian has no vanishing eigenvalues at the critical points.

Because the function $`G(\widehat{𝐮})`$ has cubic symmetry, any stationary points that occur must be members of symmetry-related families of stationary points. Thus, it is natural to focus attention on an irreducible wedge corresponding to 1/48 of the unit sphere, in terms of which Eq. (15) can be rewritten as

$$\underset{j}{}n_j\gamma _j=2$$
(16)
where $`j`$ runs over the stationary points located in the interior or on the boundary of the irreducible wedge, $`n_j`$ is a degeneracy factor counting the number of images of the stationary point generated by the cubic symmetry group operations, and $`\gamma _j`$ is $`+1`$ for a maximum or minimum and $`1`$ for a saddle point. Symmetry requires only that there be stationary points at T, O, and R; those at T and R must be maxima or minima while that at O can be of any type. Such a symmetry set of stationary points may or may not be large enough to satisfy the Morse relations as well. A set which contains the smallest number of critical points satisfying both the topological and symmetry requirements is denoted a minimal set.

This type of analysis is illustrated in Fig. 1. Figure 1(a) shows an arbitrarily chosen function $`G(\widehat{𝐮})`$ corresponding to $`c_4=c_6=0`$ and $`c_8=1`$ in Eq. (11), having maxima at T, R, and O; saddle points on the lines connecting T–O and O–R; and a minimum on the segment connecting T–R. Figure 1(b) specifies the notation that we shall use to identify the symmetry points and lines. These are chosen to correspond to the labels of distorted crystal structures (tetragonal, rhombohedral, or orthorhombic for the symmetry points; monoclinic of type ‘A,’ ‘B,’ or ‘C’ for the symmetry lines; and triclinic for the case of $`\widehat{𝐮}`$ pointing to the interior of the wedge). The corresponding degeneracy factors $`n`$ are given in Fig. 1(c). The overall behavior of the function inside the wedge is summarized in Fig. 1(d), in which the stationary points are marked by symbols that are open, shaded, or filled for maxima, saddle points, or minima, respectively. In this particular example, it can easily be seen that Eq. (16) is satisfied (6+8+12$``$24$``$24+24=2).

Using Eq. (16), one can begin enumerating the possible topological diagrams for the cubic system. Figs. 2(a) and 2(d) show the two possible diagrams in which there are only three stationary points in the irreducible wedge, the only possible minimal sets. In Fig. 2(a) the function has a unique minimum at T, so the crystal ground state would be tetragonal. A similar situation holds for the rhombohedral case of Fig. 2(d). Figures 2(b-c) and (e-f) show four of the six configurations that can exist when there are exactly four stationary points in the wedge. (Two others, similar to 2(b) and 2(e) but with the the role reversals T$``$R and MB$``$MC, are not shown.) In Fig. 2(c) there are two local minima, so that the ground state could be of type T or R depending on which minimum is deeper; a first-order transition between T and R may occur by a crossing of the minima. The simplest configuration leading to an orthorhombic ground state is shown in Fig. 2(f).

In order to arrive at configurations corresponding to monoclinic or triclinic phases, more than four stationary points are required. Figs. 3(a-c) illustrate possible configurations for monoclinic phases of type A, B, and C, in which the order parameter lies on the T–R, R–O, and O–T lines, respectively. Figure 3(d) illustrates a possible triclinic phase. It is straightforward to check that Eq. (16) is satisfied for each configuration in Figs. 2 and 3.

Of course, if the expansion of Eq. (8) or (11) is truncated at a certain order, the possible types of topological behavior will be limited by the enforced “smoothness” of the functions allowed at that order. The purpose of the following Section is to explore precisely this issue, i.e., to clarify what types of phases and phase transitions can occur at each order in the expansion. With this information in hand, one can then easily determine what is the minimal model needed to study a particular physical phenomenon of interest.

## IV Results

### A Fourth-order theory

If the expansion (11) is truncated at fourth order, the only non-constant cubic invariant is $`f_4`$ of Eq. (14). The T or R phase is favored as in Fig. 2(a) or 2(b) if $`c_4>0`$ or $`c_4<0`$, respectively. A transition between T and R phases occurs at $`c_4=0`$, but this transition is unphysical because the energy surface is perfectly flat at the transition. This degenerate behavior is an artifact of the truncation to fourth order.

### B Sixth-order theory

At sixth order in the expansion (11), the behavior is governed by the two coefficients $`c_4`$ and $`c_6`$. Clearly a common rescaling of the coefficients by a positive scale factor is irrelevant, so without loss of generality we may set $`c_4^2+c_6^2=1`$. It is then convenient to let

$$c_4=\mathrm{cos}(\alpha ),c_6=\mathrm{sin}(\alpha ).$$
(17)
By doing numerical calculations and making plots such as that shown in Fig. 1(a), we find that only three phases exist at this level. For $`\mathrm{tan}^1(3)<\alpha <\pi /2`$, the system is in the T ground state; for $`\pi /2<\alpha <5\pi /4`$, the system is in the R ground state; and for $`5\pi /4<\alpha <2\pi \mathrm{tan}^1(3)`$ the system is in the O ground state. Fig. 2 illustrates the sequence of topologies traversed as $`\alpha `$ is increased. Starting from $`\alpha =\mathrm{tan}^1(3)`$, one finds T phases corresponding to Figs. 2(a-c) consecutively, until a first-order transition occurs to the R state at $`\alpha =\pi /2`$. Then one finds R phases as illustrated in Figs. 2(c-e) until a first-order transition occurs to the O phase at $`\alpha =5\pi /4`$. Finally, O phases corresponding to Figs. 2(e-f) are found as $`\alpha `$ is increased further up to $`2\pi \mathrm{tan}^1(3)`$. The transition from the O to the T phase is degenerate in the sense that the energy surface becomes exactly flat along the entire MC symmetry line at the critical $`\alpha `$ (this being an artifact of truncation to low order).

It is important to emphasize that that no monoclinic phase is possible in the sixth-order model. In the range $`\pi +\mathrm{tan}^1(3/2)<\alpha <2\pi \mathrm{tan}^1(3)`$ (i.e., $`1.3128\pi <\alpha <1.6024\pi `$) in which such a 1D local minimum appears along the MA symmetry line, it is always unstable in the second dimension (i.e., it is a saddle point), and the true minimum is at the O point as in Fig. 2(f).

Further details of the sequence of topologies and the boundaries between them is given in the Appendix.

### C Eighth-order theory

When the model of Eq. (11) is carried to eighth order, it no longer becomes profitable to enumerate every possible topology, as was done above for the sixth-order theory. Instead, we choose to focus just on the “phase diagram” that is generated by finding the ground-state symmetry as a function of the parameters $`c_4`$, $`c_6`$, and $`c_8`$. A common scaling of the magnitudes of these coefficients is again unimportant, so we can describe the phase diagram in terms of two dimensionless parameters that we may take as

$`c_4`$ $`=`$ $`\mathrm{cos}(\alpha ),`$ (18)
$`c_6`$ $`=`$ $`\mathrm{sin}(\alpha )\mathrm{cos}(\beta ),`$ (19)
$`c_8`$ $`=`$ $`\mathrm{sin}(\alpha )\mathrm{sin}(\beta ),`$ (20)
where $`0<\alpha <\pi `$ and $`0<\beta <2\pi `$. From this point of view, the “parameter space” is just the unit sphere determined by polar and azimuthal angles $`\alpha `$ and $`\beta `$ respectively.

Figure 4 shows the phase diagram that emerges from a careful numerical study of the minimization of Eq. (11) as a function of $`\alpha `$ and $`\beta `$. The plot is a mapping of the unit sphere onto the page. The points at the “north” and “south” poles ($`\alpha =0`$ and $`\alpha =\pi `$) are the only ones accessible in the fourth-order theory (Sec. IV A); the dotted vertical lines at $`\beta =0`$ and $`\beta =\pi `$ correspond to the locus of points in parameter space that were explored by the sixth-order model (Sec. IV B).

As can be seen from Fig. 4, six of the seven possible phases (i.e., possible symmetries of a non-zero order parameter) are accessed by the eighth-order model. In addition to the T, R, and O phases that appeared already at sixth order, all three monoclinic phases (MA, MB, and MC) are now stable in some region of the phase diagram. However, the areas covered by the MC and especially MB regions are relatively small, so these phases may be harder to find in real systems than the MA phase.

Solid and dashed phase boundaries indicate transitions of first and second order, respectively, as determined numerically. According to Landau theory, there are two necessary conditions for a transition to be of second order: (i) the symmetry group of one phase must be a subset of the symmetry group of the other, $`𝒢𝒢_0`$; and (ii) from the order-parameter displacements that lead from the high-symmetry to the low-symmetry phase, it should be impossible to construct a third-order invariant of $`𝒢_0`$.  Transitions of type T–R, T–O, and O–R are necessarily first-order because of condition (i), and transitions from the R phase to MA or MB are first-order because of (ii). As can be seen from the figure, these are precisely the boundaries that were found to be of first order. All others are found to be of second order, with the exceptions of the MA–MB and MA–MC boundaries, which form a vertical line at $`\beta =3\pi /4`$ indicated by grey shading in the figure. Along this line, one finds a degenerate minimum connecting MA and MC phases as shown in Fig. 3(e) for $`0.259\pi <\alpha <0.640\pi `$, and a similar degenerate minimum connecting MA and MB phases for $`0.640\pi <\alpha <0.740\pi `$. The degenerate behavior can be traced to the fact that $`a_6=0`$ along the line $`\beta =3\pi /4`$. The triple point connecting the T, O, and MC phases at $`\alpha =\mathrm{tan}^1(1/3)=0.102\pi `$ and $`\beta =\pi `$ is also a point at which a degenerate minimum occurs, as shown in Fig. 3(f). The degenerate behaviors are artifacts of the eighth-order truncation, as will be explained more fully in Sec. IV D. Finally, the reader is reminded that because the theory is based on a single polar orientational order parameter, Eq. (6), we do not have the ability to describe transitions to or from the paraelectric C phase, nor can we describe the more complex AFE or AFD phases.

The variation of some physical variable, such as temperature, composition, or pressure, will correspond to a variation of the parameters $`\alpha `$ and $`\beta `$ of the model in a way that is not easy to predict a priori. For BaTiO3 and KNbO3, for which the observed phase transition sequence is R–O–T–C with increasing temperature, it must be the case that the system traverses a roughly vertical trajectory on the phase diagram of Fig. 4, somewhere in the range $`\pi <\beta <3\pi /2`$. (Insofar as eighth-order terms are small for these systems, the trajectory should be near $`\beta =\pi `$.) The R–O–T sequence can be visualized as traversing Figs. 2(d-f) and 2(a), in that order.

On the other hand, systems such as PZT(PbZr1-xTixO3) that exhibit a morphotropic phase boundary (MPB), i.e., an R–T transition as a function of composition $`x`$, evidently cross the first-order R–T phase boundary in the vicinity of $`\beta =0`$ with increasing $`x`$. If this trajectory passes to the right of the triple point connecting R, T, and MA phases at $`\alpha =\pi /2`$, $`\beta =\mathrm{tan}^1(1/3)=0.102\pi `$ in Fig. 4, then the phase transition sequence becomes R–MA–T, as recently observed experimentally. The narrowness of the range of MA phase, only a few percent in $`x`$, suggests that the trajectory passes rather close to the triple point. In fact, there are strong experimental indications of the possible existence of a triple point in the $`x`$$`T`$ phase diagram of PZT near $`x0.47`$ and $`T100^{}`$C.  Thus, it may be that the behavior near the triple point can be explored experimentally in the PZT system.

Using the topological analysis introduced earlier, we can now clarify the nature of the T–MA–R transition sequence near the triple point. Referring to Fig. 5, we imagine traversing a downward trajectory of increasing $`\alpha `$ at fixed $`\beta `$, slightly to the right of the triple point. Starting deep in the T phase, Fig. 2(a), a saddle point detaches from the O point and traverses the MB line toward the R point, Fig. 2(b). After it passes through the R point and emerges on the “other side” (on the MA line), we find ourselves in the situation of Fig. 5(a), the R point having been converted to a local minimum. Up to this point, the global minimum remains at T. Next, the T point converts from a local minimum to a local maximum, with the simultaneous emission of a saddle point along MC and a local minimum along MA, as shown in Fig. 5(b). This event corresponds to the second-order T–MA transition. The first-order MA–R transition then occurs by the crossing of the energies of the local minima of Fig. 5(b). Once in the R phase, the local minimum and saddle point on the MA line annihilate one another to give Fig. 5(c). Finally, deep in the R phase, the MC saddle point eventually arrives at O, giving rise to the situation of Fig. 2(d).

It should thus be emphasized that the present theory makes a definite prediction about the nature of the transitions that occur in the T–MA–R transition sequence. We can predict that, if it were possible to scan with decreasing $`x`$ at a temperature below that of the triple point, one would first find a continuous rotation of the polarization from  into the (1$`\overline{1}`$0) plane starting at a critical $`x_{c2}`$, and then a discontinuous jump to the  direction when there is a crossing of the free energies of the MA and R phases at $`x_{c1}`$. Unfortunately, the fact that the T–MA and especially the R–MA boundaries lie almost vertically in the experimental $`x`$$`T`$ plane may make it difficult to test this prediction, since $`x`$ can only be varied by preparation of multiple samples. Nevertheless, this scenario seems to be supported by the numerical simulations of Ref. .

It is important to note that a triclinic phase does not occur anywhere in the phase diagram of the eighth-order model. However, we do note the possibility of observing new monoclinic phases of types MB and MC in a region near $`\beta 0.8\pi `$. In fact, the sixth-order model (vertical dotted line at $`\beta =\pi `$) comes very close to yielding a monoclinic MC phase near the triple point (grey dot at $`\alpha `$=$`\mathrm{tan}^1(3)`$, $`\beta `$=$`\pi `$) where the MC, T, and O phases are in equilibrium. If a system such as BaTiO3 or KNbO3 could somehow be perturbed so that the variation with temperature would carry the system on a trajectory passing to the left of this triple point in Fig. 4, then a novel R–O–MA–T–C (or even R–MB–O–MA–T–C) transition sequence might be observed. However, to our knowledge, no MB or MC phase has ever been observed in a cubic perovskite system.

### D Higher-order expansions

We have seen that the eighth-order expansion still does not allow for the appearance of a triclinic equilibrium phase for any parameter values. It is natural, then, to ask at what higher order in the expansion a triclinic phase can first occur. The answer is that the expansion must be carried to twelfth order before a triclinic phase can appear.

For, suppose that a triclinic phase is the ground state. Then the minimum of $`G`$ occurs at a point $`\widehat{𝐮}_0`$ located in the interior of the irreducible wedge, as illustrated in Fig. 3(d). Letting

$`\stackrel{~}{G}(\widehat{𝐮})`$ $`=`$ $`G(\widehat{𝐮})G(\widehat{𝐮}_0),`$ (21)
$`\stackrel{~}{g}_4(\widehat{𝐮})`$ $`=`$ $`g_4(\widehat{𝐮})g_4(\widehat{𝐮}_0),`$ (22)
$`\stackrel{~}{g}_6(\widehat{𝐮})`$ $`=`$ $`g_6(\widehat{𝐮})g_6(\widehat{𝐮}_0),`$ (23)
the expansion Eq. (8) can be rewritten

$$\stackrel{~}{G}(\widehat{𝐮})=\stackrel{~}{a}_8\stackrel{~}{g}_4^2+\stackrel{~}{a}_{10}\stackrel{~}{g}_4\stackrel{~}{g}_6+\stackrel{~}{a}_{12}\stackrel{~}{g}_4^3+\stackrel{~}{a}_{12}^{}\stackrel{~}{g}_6^2+\mathrm{},$$
(24)
where the $`\stackrel{~}{a}_n`$ are trivial linear combinations of the $`a_n`$. To drop the $`\stackrel{~}{a}_4`$ and $`\stackrel{~}{a}_6`$ terms, we have used the fact that $`G(\widehat{𝐮})`$ must be stationary at $`\widehat{𝐮}_0`$. (For this we also need that the gradients of $`g_4`$ and $`g_6`$ never vanish or become parallel at an interior point of the irreducible wedge; this is straightforward to confirm.)

It is now evident that if expansion (24) is truncated at eighth order, then point $`\widehat{𝐮}_0`$ is not an isolated minimum. Instead, it belongs to a degenerate locus of minima corresponding to $`\stackrel{~}{g}_4=0`$, i.e., to a contour of the function $`g_4(\widehat{𝐮})`$. This is the situation illustrated in Fig. 3(e). It occurs when $`a_6=0`$, i.e., when $`c_6+c_8=0`$, corresponding to the grey vertical line at $`\beta =3\pi /4`$ in Fig. 4.

If the expansion is carried to tenth order, then it is clear from Eq. (24) that $`\stackrel{~}{G}`$ still vanishes on this same contour. The two-dimensional Hessian matrix $`H_{\mu \nu }=d^2G/d\widehat{u}_\mu d\widehat{u}_\nu `$ then takes the form

$$H=\left(\begin{array}{cc}0& d\\ d& e\end{array}\right),$$
(25)
where indices $`\mu `$=1 and 2 correspond to the directions parallel and perpendicular to the $`g_4`$ contour, respectively. $`d0`$ as long as $`\stackrel{~}{a}_{10}0`$, in which case $`detH<0`$. Thus, at tenth order, the stationary point $`\widehat{𝐮}_0`$ cannot be a minimum; instead, it is generally an isolated saddle point.

Finally, it is clear that the point $`\widehat{𝐮}_0`$ can be a local minimum if the expansion is carried to twelfth order. For suppose $`\stackrel{~}{a}_{10}`$=$`\stackrel{~}{a}_{12}`$=0, and $`\stackrel{~}{a}_8`$ and $`\stackrel{~}{a}_{12}^{}`$ are positive. Then $`\stackrel{~}{G}=\stackrel{~}{a}_8\stackrel{~}{g}_4^2+\stackrel{~}{a}_{12}^{}\stackrel{~}{g}_6^2`$ is obviously positive definite, and $`\widehat{𝐮}_0`$ is a true isolated (global) minimum.

Concluding this section, we find that it is necessary to go to surprisingly high order in the free-energy expansion in order to stabilize a ferroelectric state in which there are no symmetry constraints on the order parameter $`𝐏`$. Specifically, we find that cubic invariants of twelfth or higher order have to be included to to stabilize such a triclinic phase. We conclude that the discovery (or synthesis) of a material having such behavior may be challenging, but is by no means impossible.

## V A microscopic model

When structural transitions have some order-disorder character, a model free energy (e.g., (5) or (11)) expressed as a function of a macroscopic order parameter (e.g., $`𝐮`$ or $`\widehat{𝐮}`$) provides little insight into the local structural fluctuations that underlie the transitions. In such a case, a more appropriate microscopic picture of the high-symmetry phase may be one in which local regions have undergone a symmetry-lowering structural distortion, but in such a way that long-range order has not set in. For example, the description of the phase transition sequence of BaTiO3 and KNbO3 in terms of the well-known “eight-site model” assumes the presence of random local rhombohedral displacements in the orthorhombic, tetragonal, and cubic phases. In the present case, Noheda et al. have concluded from their own structural analysis of PZT and that of Corker et al. of rhombohedral PZT, that there may be random local monoclinic displacements which order variously to yield the tetrahedral, rhombohedral, or monoclinic phases near the morphotropic phase boundary. One could then describe the system in terms of fluctuations between minima of a “24-site model.”

To make these ideas more precise, suppose that the local displacements of Noheda et al. arise from an optical branch of the phonon spectrum, and let $`𝐮_l`$ be the vector “local mode amplitude” for that branch within the $`l`$th unit cell. Take as a model of the free energy

$$F[𝐮]=\underset{l}{}f(𝐮_l)+\frac{1}{2}\underset{lm}{}^{}f_{lm}(𝐮_l,𝐮_m).$$
(26)
The in-cell energy $`f`$ is presumed to be strongly nonlinear; the inter-cell coupling $`f_{lm}`$ may either be bilinear in $`𝐮_l`$ and $`𝐮_m`$, or of higher order. $`f`$ is also presumed to dominate the $`f_{lm}`$ so that, to a good approximation, its global minima establish the possible magnitude and orientations $`𝐮_l^\alpha `$ of $`𝐮_l`$. $`f(𝐮_l)`$ can be subjected to precisely the same methods of analysis as applied to $`F(𝐮)`$ in Secs. II-IV above, with parallel results. The $`𝐮_l^\alpha `$ so obtained can then be substituted into $`F`$ in Eq. (26), yielding.

$$F[\alpha ]=F_0+\frac{1}{2}\underset{lm}{}^{}f_{lm}(𝐮_l^\alpha ,𝐮_m^\beta ).$$
(27)
This can be regarded as a 24-site version of the $`q`$-state Potts model,

A statistical analysis can be carried out for various forms of $`f_{lm}`$ (e.g., bilinear) to capture the cubic (completely disordered), tetragonal (partially disordered), rhombohedral (partially disordered), and monoclinic (fully ordered) phases of the case where $`𝐮_l^\alpha `$ takes on all of the 24 symmetrically equivalent monoclinic displacements. This type of analysis has already been carried out for orientational order-disorder transitions and plastic crystals.

## VI Summary and conclusions

The original Devonshire theory gave a natural explanation for the appearance of tetragonal, orthorhombic, and rhombohedral phases in materials such as BaTiO3 based on a sixth-order free-energy expansion. Here, we have confirmed that these ferroelectric phases, in which the order parameter $`𝐏`$ is confined to a symmetry axis, are the only ones permitted by the sixth-order version of the theory. Moreover, we have clarified the nature of the phases that may be expected to appear at higher orders in the expansion. In particular, we have shown that the extension of the theory to eighth order allows one to describe, in addition, three kinds of monoclinic phase in which $`𝐏`$ is confined only to a symmetry plane. To obtain a triclinic phase in which $`𝐏`$ is unconstrained by symmetry, we have shown that a twelfth-order version of the theory is needed. A topological analysis of the critical points of the energy surface has been used to facilitate the discussion of the relevant phases and phase transitions.

The present theory may provide some added insight into the phase behavior of conventional ferroelectrics such as BaTiO3, but the principal new results concern cases in which the eighth-order terms are important. In particular, the theory provides a natural explanation for the monoclinic MA phase recently observed experimentally in PZT. It also predicts that if a triple point of equilibrium between T, R, and MA phases occurs, then it will be one at which first-order R–T and R–MA boundaries meet a second-order T–MA boundary. That is, $`𝐏`$ will rotate into the (1$`\overline{1}`$0) mirror plane continuously from the T side but discontinuously from the R side.

Noheda et al. have shown such a triple point in their Fig. 6, in agreement with our analysis. However, they also show a region of coexistence of the T and MA phases in the same phase diagram. In a homogeneous sample, such a coexistence region can be due to hysteresis arising from nucleation barriers to a first-order phase transition. The eighth-order theory predicts the T–MA phase boundary to be of second order, in which case there can be no hysteresis or nucleation barriers. In principle, higher-order contributions to the free energy could be large enough to change the order of the T–MA transition. However, a more likely explanation is that the samples studied by Noheda et al. may be inhomogeneous. We note that these authors did not report hysteresis but did report a two-phase coexistence region near the R–T boundary between 500 and 575 K in Fig. 5 of Ref.  Attributing this to inhomogeneity, we estimate that the concentration variation may be of order 1% from the slope of the R–T/M–T phase boundary in Fig. 6 of Ref. . Composition inhomogeneity of that magnitude would be sufficient to account for the T–MAphase coexistence shown in that figure.

Finally, the work may provide some guidance in the search for even more novel MC, MB, and triclinic ferroelectric phases. It also may be of utility in other kinds of cubic systems with other kinds of vector order parameters, e.g., ferromagnetic systems.

###### Acknowledgements.

D.V. acknowledges support of ONR Grant N00014-97-1-0048. We wish to thank L. Bellaiche, A. García, A. Khachaturyan, and B. Noheda for useful comments and discussions.

## A Details of sixth-order theory

The purpose of this Appendix is to give further details about the sequence of transitions that occurs in the sixth-order model of Sec. IV B. Recall that the behavior in this model is governed by a single dimensionless parameter $`\alpha `$ defined via Eq. (17).

There are ten critical values of $`\alpha `$ that we can define as

$`\alpha _1`$ $`=`$ $`0.3976\pi =\mathrm{tan}^1(3)`$ (A1)
$`\alpha _2`$ $`=`$ $`0.2256\pi =\mathrm{tan}^1(6/7)`$ (A2)
$`\alpha _3`$ $`=`$ $`0.25\pi `$ (A3)
$`\alpha _4`$ $`=`$ $`0.3128\pi =\mathrm{tan}^1(3/2)`$ (A4)
$`\alpha _5`$ $`=`$ $`0.5\pi `$ (A5)
$`\alpha _6`$ $`=`$ $`0.6024\pi =\pi \mathrm{tan}^1(3)`$ (A6)
$`\alpha _7`$ $`=`$ $`1.2256\pi =\pi +\mathrm{tan}^1(6/7)`$ (A7)
$`\alpha _8`$ $`=`$ $`1.25\pi `$ (A8)
$`\alpha _9`$ $`=`$ $`1.3128\pi =\pi +\mathrm{tan}^1(3/2)`$ (A9)
$`\alpha _{10}`$ $`=`$ $`1.5\pi `$ (A10)
The sequence of phases can be followed on Fig. 4 by tracing the vertical dotted lines, first from top to bottom at $`\beta =0`$ for $`0<\alpha <\pi `$, and then from bottom to top at $`\beta =\pi `$ for $`\pi <\alpha <2\pi `$.

In the T phase, the system exhibits the topology of Fig. 2(a) for $`\alpha _1<\alpha <\alpha _2`$; Fig. 2(b) for $`\alpha _2<\alpha <\alpha _4`$; and Fig. 2(c) for $`\alpha _4<\alpha <\alpha _5`$. At $`\alpha _2`$, point O converts to a local maximum and simultaneously a saddle point appears along the MB symmetry line. There is an irrelevant crossing of the R and O maxima at $`\alpha _3`$. At $`\alpha _4`$, R becomes a local minimum, and the saddle point switches from the MB to the MA symmetry line. The transition from T to R is first order at $`\alpha _5`$.

The system falls into an R ground state corresponding to Fig. 2(c) for $`\alpha _5<\alpha <\alpha _6`$; Fig. 2(d) for $`\alpha _6<\alpha <\alpha _7`$; and Fig. 2(e) for $`\alpha _7<\alpha <\alpha _8`$. The MA saddle point disappears and T is converted to a maximum at $`\alpha _6`$, and a new MB saddle point emerges with the conversion of O to a local minimum at $`\alpha _7`$. The transition from R to O is first order at $`\alpha _8`$.

Finally, an O phases occurs as illustrated in Fig. 2(e) for $`\alpha _8<\alpha <\alpha _9`$, and as in Fig. 2(f) for $`\alpha _9<\alpha <\alpha _1+2\pi `$. The saddle point at MB vanishes and R is converted to a local maximum at $`\alpha _9`$. There is an irrelevant crossing of the T and R maxima at $`\alpha _{10}`$. The transition from O to T at $`1.6024\pi `$ is singular, in that the energy surface becomes exactly flat along the entire MC symmetry line.


format: text

text:
# Franklin’s argument proves an identity of Zagier

## 1 Introduction

We use the standard $`q`$-series notation:

$$(a)_n=\underset{k=1}{\overset{n}{}}(1aq^{k1})$$

where $`n`$ is a nonnegative integer or $`n=\mathrm{}`$. Euler’s pentagonal number theorem states that

$$(q)_{\mathrm{}}=1+\underset{r=1}{\overset{\mathrm{}}{}}(1)^r(q^{r(3r1)/2}+q^{r(3r+1)/2}).$$
(1)
Recently Zagier proved the following remarkable identity

###### Theorem 1

$$\underset{n=0}{\overset{\mathrm{}}{}}[(q)_{\mathrm{}}(q)_n]=(q)_{\mathrm{}}\underset{k=1}{\overset{\mathrm{}}{}}\frac{q^k}{1q^k}+\underset{r=1}{\overset{\mathrm{}}{}}(1)^r[(3r1)q^{r(3r1)/2}+3rq^{r(3r+1)/2}].$$
(2)
This is \[4, Theorem 2\] slightly rephrased.

Equation (1) has a combinatorial interpretation. The coefficient of $`q^N`$ in $`(q)_{\mathrm{}}`$ equals $`d_e(N)d_o(N)`$ where $`d_e(N)`$ (respectively $`d_o(N)`$) is the number of partitions of $`N`$ into an even (respectively odd) number of distinct parts. Franklin  showed that

$$d_e(N)d_o(N)=\{\begin{array}{cc}(1)^r& \text{if }N=\frac{1}{2}r(3r\pm 1)\text{ for a positive integer }r\text{,}\hfill \\ 0& \text{otherwise.}\hfill \end{array}$$

His proof was combinatorial. He set up what was almost an involution on the set of partitions of $`N`$ into distinct parts. This “involution” reverses the parity of the number of parts. However there are certain partitions for which his map is not defined. These exceptional partitions occur precisely when $`N=\frac{1}{2}r(3r\pm 1)`$, and so account for the nonzero terms on the right of (1).

We show that Zagier’s identity has a similar combinatorial interpretation, which, miraculously, Franklin’s argument proves at once.

## 2 Proof of Theorem 1

We begin by recalling Franklin’s “involution”. Let $`𝒟_N`$ denote the set of partitions of $`N`$ into distinct parts and let $`𝒟=_{N=0}^{\mathrm{}}𝒟_N`$. For $`\lambda 𝒟_N`$ let $`N_\lambda =N`$, $`n_\lambda `$ be the number of parts in $`\lambda `$ and $`m_\lambda `$ be the largest part of $`\lambda `$ (if $`\lambda `$ is the empty partition of 0 let $`m_\lambda =0`$). Then

$$(q)_{\mathrm{}}=\underset{\lambda 𝒟}{}(1)^{n_\lambda }q^{N_\lambda }.$$
(3)
Let $`\lambda `$ be a non-empty partition in $`𝒟`$. Denote its smallest part by $`a_\lambda `$. If the parts of $`\lambda `$ are $`\lambda _1>\lambda _2>\lambda _3>\mathrm{}`$ let $`b=b_\lambda `$ denote the largest $`b`$ such that $`\lambda _b=\lambda _1+1b`$ (so that $`\lambda _k=\lambda _1+1k`$ if and only if $`1kb`$). If $`\lambda 𝒟`$ is not *exceptional* (we shall explain this term shortly), then we define a new partition $`\lambda ^{}`$ as follows. If $`a_\lambda b_\lambda `$ we obtain $`\lambda ^{}`$ by removing the smallest part from $`\lambda `$ and then adding 1 to the largest $`a_\lambda `$ parts of this new partition. If $`a_\lambda >b_\lambda `$ we obtain $`\lambda ^{}`$ by subtracting 1 from the $`b_\lambda `$ largest parts of $`\lambda `$ and then appending a new part $`b_\lambda `$ to this new partition.

The exceptional partitions are those for which this procedure breaks down. We regard the empty partition as exceptional, also we regard those for which $`n_\lambda =b_\lambda `$ and $`a_\lambda =b_\lambda `$ or $`b_\lambda +1`$. If $`\lambda `$ is not exceptional, then neither is $`\lambda ^{}`$ and $`\lambda ^{\prime \prime }=\lambda `$ and $`(1)^{n_\lambda ^{}}=(1)^{n_\lambda }`$. Thus on the right side of (3) the contributions from non-exceptional partitions cancel. The non-empty exceptional partitions are of two forms: for each positive integer $`r`$ we have $`\lambda =(2r1,2r2,\mathrm{},r+1,r)`$ for which $`n_\lambda =r`$, $`m_\lambda =2r1`$ and $`N_\lambda =\frac{1}{2}r(3r1)`$, and we have $`\lambda =(2r,2r1,\mathrm{},r+2,r+1)`$ for which $`n_\lambda =r`$, $`m_\lambda =2r`$ and $`N_\lambda =\frac{1}{2}r(3r+1)`$. Thus from (3) we deduce (1).

If $`\lambda 𝒟`$ is non-exceptional, then either $`n_\lambda ^{}=n_\lambda 1`$ in which case $`m_\lambda ^{}=m_\lambda +1`$ or $`n_\lambda =n_\lambda +1`$ in which case $`m_\lambda ^{}=m_\lambda 1`$. In each case $`m_\lambda ^{}+n_\lambda ^{}=m_\lambda +n_\lambda `$. It follows that in the sum

$$\underset{\lambda 𝒟}{}(1)^{n_\lambda }(m_\lambda +n_\lambda )q^{N_\lambda }$$

the terms corresponding to non-exceptional $`\lambda `$ cancel and so we get only the contribution from exceptional $`\lambda `$. Thus

$$\underset{\lambda 𝒟}{}(1)^{n_\lambda }(m_\lambda +n_\lambda )q^{N_\lambda }=\underset{r=1}{\overset{\mathrm{}}{}}(1)^r[(3r1)q^{r(3r1)/2}+3rq^{r(3r+1)/2}].$$
(4)
This sum occurs in (2), which will follow by analysing the left side of (4).

We break this into two sums. The first

$$\underset{\lambda 𝒟}{}(1)^{n_\lambda }m_\lambda q^{N_\lambda }$$

is dealt with in \[2, Theorem 5.2\]. We repeat their argument. The coefficient of $`q^N`$ in $`(q)_{\mathrm{}}(q)_n`$ is the sum of $`(1)^{n_\lambda }`$ over all $`\lambda 𝒟_N`$ having a part strictly greater than $`n`$. Such a $`\lambda `$ is counted for exactly $`m_\lambda `$ different $`n`$ so that

$$\underset{n=0}{\overset{\mathrm{}}{}}[(q)_{\mathrm{}}(q)_n]=\underset{\lambda 𝒟}{}(1)^{n_\lambda }m_\lambda q^{N_\lambda }.$$
(5)
For each positive integer $`k`$,

$$\frac{q^k}{1q^k}(q)_{\mathrm{}}=(1q)(1q^2)\mathrm{}(1q^{k1})(q^k)(1q^{k+1})\mathrm{}.$$

The coefficient of $`q^N`$ in this product is the sum of $`(1)^{n_\lambda }`$ over all $`\lambda 𝒟_N`$ having $`k`$ as a part. Such a $`\lambda `$ occurs for $`n_\lambda `$ distinct $`k`$, and summing we conclude that

$$(q)_{\mathrm{}}\underset{k=1}{\overset{\mathrm{}}{}}\frac{q^k}{1q^k}=\underset{\lambda 𝒟}{}(1)^{n_\lambda }n_\lambda q^{N_\lambda }.$$
(6)
Combining (4), (5) and (6) gives (2).

## 3 Another identity

We can prove another identity by the same method. As before Franklin’s argument proves that

$$\underset{\lambda 𝒟}{}(1)^{n_\lambda }x^{m_\lambda +n_\lambda }q^{N_\lambda }=1+\underset{r=1}{\overset{\mathrm{}}{}}(1)^r[x^{3r1}q^{r(3r1)/2}+x^{3r}q^{r(3r+1)/2}].$$
(7)
We now give an alternative way of expressing the left side of (7). Consider the contribution due to the partitions $`\lambda `$ with $`m_\lambda `$ equalling a fixed $`m1`$. The product

$$(xq)_m=\underset{j=1}{\overset{m}{}}(1xq^j)=\underset{\genfrac{}{}{0pt}{}{\lambda 𝒟}{m_\lambda m}}{}(1)^{n_\lambda }x^{n_\lambda }q^{N_\lambda }$$

and so

$$\underset{\genfrac{}{}{0pt}{}{\lambda 𝒟}{m_\lambda =m}}{}(1)^{n_\lambda }x^{m_\lambda +n_\lambda }q^{N_\lambda }=x^m[(xq)_m(xq)_{m1}].$$

Hence

$`{\displaystyle \underset{\lambda 𝒟}{}}(1)^{n_\lambda }x^{m_\lambda +n_\lambda }q^{N_\lambda }`$ $`=`$ $`1+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}x^m[(xq)_m(xq)_{m1}]`$
$`=`$ $`(1x){\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}(xq)_mx^m`$
$`=`$ $`{\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}(x)_{r+1}x^r.`$
Hence

$$\underset{r=0}{\overset{\mathrm{}}{}}(x)_{r+1}x^r=1+\underset{r=1}{\overset{\mathrm{}}{}}(1)^r[x^{3r1}q^{r(3r1)/2}+x^{3r}q^{r(3r+1)/2}].$$
(8)
This can also be proved by noting that both sides satisfy the recurrence $`S(x)=1qx^2s^2x^3S(qx)`$. This appears as exercise 10 in Chapter 2 of  whose solution is outlined in . Zagier deduces (2) from (8), essentially by carefully differentiating with respect to $`x`$ and setting $`x=1`$.

## 4 Acknowledgments

The author wishes to thank George Andrews and Don Zagier for supplying him with copies of  and .


text:
# Intermediate Velocity Clouds – Distances, Infall and Shadowing

## 1. Distances and Infall

For each IVC, a set of target stars was observed spectroscopically in the Na i doublet region for evidence of absorption from the target IVC (details in Gladders et al. 1998). The distances and derived cloud masses are presented in Table 1.

The IVC distances permit us to compare cloud velocities with Galactic infall models. The general trend of the estimated infall velocities appears consistent with an increase in velocity with increasing height; this trend is not expected for clouds following ballistic trajectories (see Clarke et al. 1999 for details).

## 2. Shadowing

The nature of the SXRB is in dispute. One model suggests that most of the SXRB arises locally from the low density plasma surrounding the Sun and extending for a few hundred parsecs (Sanders et al. 1977). Other models propose that most of the SXRB originates at large distances in a Galactic halo (e.g. Marshall and Clark 1984).

Soft X-ray background (SXRB) measurements toward Draco show a reduction of $``$ 50 % in the 1/4 keV emission relative to the adjacent sky, consistent with photoelectric absorption of the SXRB by Draco (Burrows and Mendenhall 1991; Snowden et al. 1991). Furthermore, the profile of the cloud’s dust emission (Figure 1a) is closely followed by the X-ray shadowing (Figure 1b). Our distance to Draco places it beyond the local plasma and therefore the X-ray emission it is shadowing must have an origin at larger Galactic distances. Specifically, 50 % or more of the SXRB emission originates from distances greater than 460 pc. Note that a sample of X-ray shadowing clouds at known distances (such as the DDO IVC sample) can be used to infer the detailed distribution of the SXRB through determining the foreground and background components (work in progress).

### Acknowledgments.

We acknowledge the work of the entire DDO IVC Distance Project team and thank the organizers of this meeting.

Refernces

Burrows, D.N., & Mendenhall, J.A., Nature, 351, 629
Clarke, T.E., et al. 1999, ASP Conference Series 168
Gladders, M.D. et al., 1998, ApJ, 507, 161
Marshall, F.J., & Clark, G.W. 1984, ApJ, 287, 633
Sanders, W.T., et al. 1977, ApJL, 217, 87
Snowden, S.L., et al. 1991, Science, 525, 1529


text:
# Reprocessed emission line profiles from dense clouds in geometrically thick accretion engines

## 1 Introduction

Active galactic nuclei (AGN) and some X-ray binaries are thought to be powered by accretion onto massive black holes (e.g. Rees 1984; Tanaka 2000). Observations suggest that cold matter in the accretion engine reprocesses the primary X-ray continuum (e.g. George and Fabian 1991). In particular, a fluorescent iron K$`\alpha `$ emission line profile is imprinted on the spectrum near 6.4 keV. The iron line profile is shaped by Doppler and gravitational frequency shifts (Fabian et al. 1989), and hence contains information about geometry and dynamics in the innermost regions of the accretion engine (see Fabian et al. 2000 for a comprehensive review).

Most modelling of iron line profiles has been carried out within the paradigm of geometrically thin, optically thick discs. The initial flat disc model of Fabian et al. (1989) was generalised to Kerr black holes by Laor (1991). More recently, Pariev and Bromley (1998) have considered effects of a small disc thickness and turbulence.

Thin disc line profiles successfully reproduce many features in the experimental data (e.g. Tanaka et al. 1995; Nandra et al. 1997). However there are some observations which remain difficult to understand within this paradigm. For example, a preference for low inclinations, even amongst Seyfert 2s (Turner et al. 1998), disagreements between the predicted inclination from iron line modelling and the predicted inclination from other phenomena (Nishiura et al. 1998; Sulentic et al. 1998a; Wang et al. 1999), and problems with the observed reprocessed fraction being too low or not correlated with line width (e.g. Lee et al. 2000; Chiang et al. 2000; Weaver 2000; Sulentic et al. 1998b).

Such puzzles may be resolved within the paradigm of thin discs, for instance by introducing additional emission components from an obscuring torus (Weaver and Reynolds 1998) or invoking flared or warped discs in which the inner and outer regions are not aligned (Blackman 1999; Hartnoll and Blackman 2000). In any case, alternative plausible geometries still merit exploration to provide proper diagnostics for comparison.

A particularly interesting alternative extends from work on the survival of cold clouds capable of reprocessing X-ray flux in the AGN engines (e.g. Ferland and Rees 1988; Guilbert and Rees 1988; Celotti, Fabian and Rees 1992). Physical constraints on the size of the dense clouds have been studied by Kuncic et al. (1996). They found a sizeable parameter regime in which cool, dense, small clouds can coexist with a hot medium in a dynamical equilibrium. Though X-ray continuum features of the embedded cloud model have been studied by several authors (Sivron and Tsuruta 1993; Bond and Matsuoka 1993; Nandra and George 1994; Collin-Souffrin et al. 1996), detailed calculations of reprocessed line profiles have not been carried out. Karas et al. (2000) have studied line emission from reprocessing by a nearly spherical shell of clouds at a specific radius that results from a standard thin disc being disrupted due to instabilities. They give resulting line profiles from illumination by a central source. This is a different scenario from what we will consider here. Other related calculations of note include Cao and Zhang (1991), who calculated line profiles for radially infalling gas.

In this paper we calculate the (iron) line profiles that result from various distributions of small, cold clouds embedded in an optically thin, geometrically thick accretion disc or corona. As well as illumination from a central source, we include an illumination by an ‘X-ray bath’ of radiation; each cloud is taken to be small enough to be illuminated isotropically at their positional location by this bath component. This two component illumination is motivated from considering geometrical thick models such as, for example, advection dominated accretion flows (ADAFs) (e.g. Rees et al. 1982; Narayan and Yi 1994; see Narayan et al. 1998 for a review), which can be interpreted to have a synchrotron + Bremsstrahlung centrally dominated illumination (Mahadevan 1997) and a Bremsstrahlung X-ray bath. Depending on the accretion rate of the disc, the synchrotron or Bremsstrahlung will dominate. ADAFs predict luminosities significantly less than the Eddington value $`L_E`$, and may just be applicable to Seyferts, with typical luminosities $`L_E/100`$.

We will use several flow solutions to prescribe a velocity field for the clouds; standard full ADAF solutions, pure radial infall, and bipolar outflows. The profiles calculated may also therefore be relevant for the dynamics of convection dominated accretion flows (CDAFs) (Quataert and Gruzinov 2000; Narayan et al. 2000), Advection Dominated Inflow Outflow Models (ADIOS) (Blandford and Begelman 1999) or other thick disc models. The microphysics of the direct X-ray line emission for ADAFs have been studied in non-relativistic flows (e.g. Narayan and Raymond 1999; Perna et al. 2000), but not the presence of cold reprocessing matter, nor the associated relativistic iron flourescence line profiles. Like Karas et al. (2000), we take the clouds to be optically thick, and allow the clouds to shadow emission from other clouds.

Some profiles for relativistic outflows have been calculated by Wang et al. (2000). These outflows are different from those considered here in that they are pencil beams along the spin axis and not radial outflow in a geometrically thick disc. The spherical geometry, and covering by clouds, can avoid the large blue shifted component that has not been observed in iron lines and that is characteristic of outflows.

It has been suggested (e.g. Nandra 2000) that a geometrically thick distribution of clouds would result in more blue-shifted emission than is observed. The data from a variety of sources are also not inconsistent with a somewhat larger blue wing (Nandra et al. 1997). But we find that because models of geometrically thick discs can generally have smaller bulk velocities in the central regions than the flat disc keplerian velocities, smaller Doppler blue shifts result, and gravitational redshift becomes the dominant spectral feature. For these reasons also, the line profiles for thick discs will in general be less inclination dependent than those from thin discs. This will be discussed in more detail below. It may help explain problems encountered with inclinations in thin disc modelling. When the velocities are faster, the profiles are very similar to the standard profiles obtained from thin disc modelling.

In section 2 we describe the disc model of dense clouds in geometrically thick discs. Section 3 contains the calculation of iron line profiles from this model. Section 4 presents the resulting line profiles and a discussion of the main features. Section 5 is the conclusion.

## 2 The embedded cloud model: dense clouds and thick discs

The model consists of small optically thick clouds embedded in, and hence following, the fluid flow in an otherwise optically thin disc. Although the clouds may have a short lifetime, the continuous creation and destruction of clouds allows the definition of an effective number density of clouds at each point, $`n_c(𝐫)`$. We consider these clouds to be both bathed in an locally isotropic radiative flux (e.g. Bremsstrahlung) and to be illuminated by a central source (e.g. synchrotron dominated). The clouds are assumed to be small compared to disc length scales (see Celotti et al. 1992 and Kuncic et al. 1996 for partial justification).

Calculating the profile then involves an optical depth calculation; we integrate the number density times the cross section along the photon trajectories to obtain a probability of obscuration. Contributions are integrated across the whole disc to obtain a profile, as is done for the usual thin disc, except that now we are integrating over three dimensions.

We consider several different solutions to the accretion flow involving optically thin, geometrically thick engines: self-similar ADAFs (Narayan and Yi 1995; Perna et al. 2000), bipolar outflows (Xu and Chen 1997; Perna et al. 2000), and radial infall (appendix B of Narayan and Yi 1995, similar to Bondi 1952).

### 2.1 self-similar ADAF

Self-similar accretion solutions are described in Narayan and Yi (1995). The approximate analytic solutions, given in Perna et al. (2000), for the density ($`\rho `$), velocities ($`v_r,v_\theta ,v_\varphi `$), and sound speed ($`c_s`$) are

$$\rho =\rho _0r^{\frac{3}{2}}[c_s^2(\theta )]^{2.25},$$
(1)
$$v_r=\frac{3\alpha }{5+2ϵ}r\mathrm{\Omega }_K(r)\mathrm{sin}^2\theta ,$$
(2)
$$v_\theta =0,$$
(3)
$$v_\varphi =\left(\frac{2ϵ}{5+2ϵ}\right)^{\frac{1}{2}}r\mathrm{\Omega }_K(r)\mathrm{sin}\theta ,$$
(4)
$$c_s=r\mathrm{\Omega }_K(r)\left[\frac{2}{5}\left(\frac{2}{5}\frac{2}{5+2ϵ}\right)\mathrm{sin}^2\theta \right]^{\frac{1}{2}},$$
(5)
where $`\mathrm{\Omega }_K(r)`$ is the Keplerian angular velocity, $`\mathrm{\Omega }_K(r)=r^{\frac{3}{2}}`$ in our units and we are using standard spherical polar coordinates, $`(r,\theta ,\varphi )`$.

There are two free parameters, the Sakura-Sunyaev viscosity parameter, $`\alpha `$, and a thermodynamic parameter,

$$ϵ=\left(\frac{5/3\gamma }{\gamma 1}\right),$$
(6)
with $`\gamma `$ the ratio of specific heats. We will consider $`\alpha =0.35`$ and $`ϵ=10,0.1`$. Solutions are taken to be in the fully advection dominated case. Note that:

1. Lower $`ϵ`$ gives lower angular speeds.
2. For $`ϵ=10`$ most of the matter is close to the equatorial plane, whilst for $`ϵ=0.1`$ the matter is essentially uniformly distributed in $`\theta `$.
3. For larger values of $`\alpha `$, such the $`\alpha =0.35`$ we consider, $`v_r`$ is very similar to $`v_\varphi `$. For smaller $`\alpha `$, the radial motion becomes less important.
4. $`c_s(\theta )`$ is smaller on the equatorial plane and is smaller for larger $`ϵ`$. This will mean that, at fixed radius, the clouds are larger on the equatorial plane and larger in high $`ϵ`$ cases.

### 2.2 Radial infall

The self-similar forms of the previous section admit a radially infalling solution (Narayan and Yi 1995, it is essentially Bondi accretion with viscosity.) We take the radial infall to be at the Keplerian velocity, $`r\mathrm{\Omega }_K(r)`$, which implies that the sound speed is $`c_s=\frac{r\mathrm{\Omega }_K(r)}{\sqrt{5}}`$. This solution corresponds to a small $`ϵ`$, and so represents a geometrically thick, optically thin, advection dominated flow. The density is taken to be constant in $`\theta `$, so $`\rho =\rho _0r^{\frac{3}{2}}`$.

### 2.3 Bipolar outflow

Outflow solutions have been considered by several authors (e.g. Xu and Chen 1997; Quartaert and Narayan 1999; Blandford and Begelman 1999). Perna et al. (2000) give approximate analytic expressions. We use their solutions with density independent of angle and with winds of $`p=\frac{1}{2}`$, where $`p`$ is defined so that $`\dot{M}r^p`$. The equations then become

$$\rho =\rho _0r^{\frac{3}{2}},$$
(7)
$$v_r=\frac{3\alpha }{5+2ϵ}r\mathrm{\Omega }_K(r)[\mathrm{cos}2\theta +\mathrm{cos}^2\theta ],$$
(8)
$$v_\theta =\frac{3\alpha }{5+2ϵ}\frac{r\mathrm{\Omega }_K(r)}{4}\mathrm{sin}2\theta ,$$
(9)
$$v_\varphi =\left(\frac{2ϵ}{5+2ϵ}\right)^{\frac{1}{2}}r\mathrm{\Omega }_K(r)\mathrm{sin}\theta ,$$
(10)
$$c_s=\frac{r\mathrm{\Omega }_K(r)}{\sqrt{5+2ϵ}}\sqrt{2ϵ\mathrm{cos}^2\theta }.$$
(11)
## 3 Calculation of the line profile

### 3.1 Total flux

The total flux observed at frequency $`\nu `$ is given by

$$F(\nu )=\delta (\nu \nu _0(𝐫))f(𝐫,\nu )d^3r,$$
(12)
where $`\nu _0(𝐫)`$ is the observed frequency of flux emitted from $`𝐫`$, calculated below, and $`I(𝐫,\nu )`$ is given by

$$f(𝐫,\nu )=\left[\frac{\nu }{\nu _e}\right]^3f_{out}(𝐫,𝐜(𝐫))V(𝐫),$$
(13)
where $`\left[\frac{\nu }{\nu _e}\right]^3`$ is the relativistic correction to the flux, with the frequency of emission, $`\nu _e=6.4\mathrm{keV}`$ for the iron line. $`f_{out}(𝐫,𝐜(𝐫))`$ is the flux reprocessed and re-emitted by clouds at $`𝐫`$ in the direction $`𝐜(𝐫)`$, defined such that the photons reach the observer (calculated below). Finally, $`V(𝐫)`$ is the visibility function which gives the probability of flux from $`𝐫`$ reaching the observer without being intercepted by an optically thick cloud. It is given by

$$V(𝐫)=\mathrm{max}(1P(𝐫),0),$$
(14)
with the probability of interception

$$P(𝐫_\mathrm{𝟎})=_a^bn_c(𝐫(\lambda ))A(𝐫(\lambda ),\frac{d𝐫}{d\lambda })\left|\frac{d𝐫}{d\lambda }\right|𝑑\lambda ,$$
(15)
where $`𝐫(a)=𝐫_\mathrm{𝟎}`$, the position of emission, and $`𝐫(b)=𝐫_{\mathrm{𝐨𝐛𝐬}}`$, the position of the observer, and $`\lambda `$ parametrises the path of the photon. $`n_c(𝐫)`$ is the number density of clouds in the disc and $`A(𝐫,\frac{d𝐫}{d\lambda })`$ is the cross sectional area of the clouds at $`𝐫`$, as seen by a photon travelling in direction $`\frac{d𝐫}{d\lambda }`$. These quantities are discussed below. The visibility function must also include the possibility of photons falling into the event horizon.

The dimensions of the disc are taken to be $`6r300`$, $`0\theta \pi `$ and $`0\varphi <2\pi `$, although there will not necessarily be a significant number density of clouds throughout this whole region (because the number density of particles in the disc itself gets smaller for larger $`r`$ and, depending on the thermodynamic parameter $`ϵ`$, gets smaller further form the equatorial plane). Also, the volume element has a general relativistic correction, so

$$d^3r=\frac{dr}{\sqrt{1\frac{2}{r}}}r^2\mathrm{sin}\theta d\theta d\varphi .$$
(16)
### 3.2 Photon paths

We use paths in the Schwarzchild metric, approximated by a first order perturbation. We take units with $`c=G=1`$ and take the mass of the central object to be $`M=1`$. Note that in this section $`𝐫`$ is the position of the photon travelling from the cloud to the observer, and $`𝐫_\mathrm{𝟎}`$ is the position of the emitting cloud in the disc. We write $`r_0=𝐫_\mathrm{𝟎}`$ and $`r=𝐫`$.

First order perturbation of the Schwarzchild null geodesic equations gives the photon path to be

$$\frac{1}{r(\mathrm{\Phi })}=\frac{\mathrm{sin}\mathrm{\Phi }}{2R}+\frac{1+C\mathrm{cos}\mathrm{\Phi }+\mathrm{cos}^2\mathrm{\Phi }}{4R^2}+𝒪(\frac{1}{R^3}).$$
(17)
Where $`R`$ and $`C`$ are constants, $`r(\mathrm{\Phi })`$ is the radial distance from the origin at the central object and $`\mathrm{\Phi }`$ is a polar angle on the plane of the trajectory. The trajectory plane is defined by the origin, the direction to the observer, $`\widehat{𝐧}`$, and the initial position vector $`𝐫_\mathrm{𝟎}`$.

The boundary conditions are that the photon starts at $`(r_0,\mathrm{\Phi }_0)`$ and ends at $`(\mathrm{},\mathrm{\Phi }_{obs})`$, written in $`(r,\mathrm{\Phi })`$ coordinates on the trajectory plane. The boundary conditions are simplified if we choose the axis on the trajectory plane that defines $`\mathrm{\Phi }`$ to be in the direction towards the observer, $`\widehat{𝐧}`$, so that $`\mathrm{\Phi }_{obs}=0`$. The constants are then given by

$$C=2,$$
(18)
$$R^\pm =\frac{r_0\mathrm{sin}\mathrm{\Phi }_0}{4}\pm \frac{r_0}{4}\sqrt{\mathrm{sin}^2\mathrm{\Phi }_0+\frac{4}{r_0}(1\mathrm{cos}\mathrm{\Phi }_0)^2}.$$
(19)
Note that the minus sign case results in $`R`$ negative, which is not physical, so we take the positive solution. Selecting this single path is effectively not considering orbits that loop around the central hole one or more times, consistent with the perturbative approach.

With this choice of axis, $`\mathrm{\Phi }_0`$ is given by

$$\mathrm{cos}\mathrm{\Phi }_0=\frac{𝐫_\mathrm{𝟎}\widehat{𝐧}}{r_0}.$$
(20)
In principle, there could be a problem inverting this expression due to ambiguities in the definition of the inverse cosine. However, these are avoided if we take the other axis on the trajectory plane to be $`\frac{𝐫_\mathrm{𝟎}(𝐫_\mathrm{𝟎}\widehat{𝐧})\widehat{𝐧}}{r_0}`$. This vector, orthogonal to $`\widehat{𝐧}`$, is such that $`𝐫_\mathrm{𝟎}`$ has positive component in this direction. Therefore, we may take $`\mathrm{\Phi }_0[0,\pi ]`$ with no ambiguities and the trajectory is given by

$$𝐫(\mathrm{\Phi })=r(\mathrm{\Phi })\left[\mathrm{cos}\mathrm{\Phi }\widehat{𝐧}+\mathrm{sin}\mathrm{\Phi }\frac{𝐫_\mathrm{𝟎}(𝐫_\mathrm{𝟎}\widehat{𝐧})\widehat{𝐧}}{r_0}\right].$$
(21)
This gives the probability of interception as

$$P(𝐫_\mathrm{𝟎})=_o^{cos^1\frac{𝐫_\mathrm{𝟎}\widehat{𝐧}}{r_0}}n(𝐫(\mathrm{\Phi }))A(𝐫(\mathrm{\Phi }),\frac{d𝐫}{d\mathrm{\Phi }})\left|\frac{d𝐫}{d\mathrm{\Phi }}\right|𝑑\mathrm{\Phi },$$
(22)
and the angle of emission is given by

$$𝐜(𝐫_\mathrm{𝟎})=\frac{d𝐫}{d\mathrm{\Phi }}(\mathrm{\Phi }=\mathrm{\Phi }_0).$$
(23)
The minus sign in the last two equations results because the photon travels in the direction of decreasing $`\mathrm{\Phi }`$.

### 3.3 Frequency shifts

The observed frequency is given in terms of the emitted frequency by

$$\nu _0=\nu _e\frac{(p_\mu U^\mu )_{obs}}{(p_\mu U^\mu )_{em}},$$
(24)
where $`p^\mu `$ is the photon 4-momentum and $`U^\mu `$ is the 4-velocity of the emitter or observer. We employ general relativity in the Schwarzchild metric, $`g_{\mu \nu }`$. Write $`p^\mu =(p^t,𝐩)^\mu `$, and $`U^\mu =\frac{dt}{d\tau }(1,𝐮)^\mu `$. The time component of photon 4-momentum in the Schwarzchild metric is

$$p^t=\frac{E}{1\frac{2}{r}},$$
(25)
Where E is constant along the photon trajectory. So requiring the photon 4-momentum to be null, $`g_{\mu \nu }p^\mu p^\nu =0`$, gives $`𝐩^2=g_{ij}p^ip^j=\frac{E}{1\frac{2}{r}}`$ to be the appropriate normalisation for the spatial components. In the Schwarzchild metric, we have

$$\frac{dt}{d\tau }=\frac{1}{\sqrt{1\frac{2}{r}\frac{u_r^2}{1\frac{2}{r}}u_\theta ^2u_\varphi ^2}}.$$
(26)
By definition,

$$p_\mu U^\mu =g_{\mu \nu }p^\mu U^\nu ,$$
(27)
where the observer has $`𝐮_{obs}=0`$ and $`r_{obs}\mathrm{}`$. The emitter has $`𝐩_{em}𝐜`$ and $`r_{em}=r_0`$. Putting all this into (24) gives

$$\nu _0(𝐫_\mathrm{𝟎})=\frac{\nu _e\sqrt{1\frac{2}{r_0}\frac{u_r^2(𝐫_\mathrm{𝟎})}{1\frac{2}{r_0}}u_\theta ^2(𝐫_\mathrm{𝟎})u_\varphi ^2(𝐫_\mathrm{𝟎})}}{1\frac{c_r(𝐫_\mathrm{𝟎})u_r(𝐫_\mathrm{𝟎})}{1\frac{2}{r_0}}c_\theta (𝐫_\mathrm{𝟎})u_\theta (𝐫_\mathrm{𝟎})c_\varphi (𝐫_\mathrm{𝟎})u_\varphi (𝐫_\mathrm{𝟎})}.$$
(28)
Here $`𝐮(𝐫_\mathrm{𝟎})`$ is the velocity of the clouds at $`𝐫_\mathrm{𝟎}`$. We take the velocity of clouds to trace the velocity of the fluid in the disc. This is given by one of the disc motions described in the previous section.

### 3.4 Distribution and sizes of clouds

For simplicity, we will assume that all of the clouds are spherical. The cross section is then independent of a photon’s direction of approach. If we consider clouds with the same total number of particles, this allows us to to write

$$A(𝐫)V^{\frac{2}{3}}(𝐫)\frac{1}{n_{p,c}^{\frac{2}{3}}(𝐫)},$$
(29)
where $`n_{p,c}`$ is the number density of particles in the cloud.

Thermal equilibrium gives a relation between the number density of particles in the cloud and the number density of particles in the disc, $`n`$,

$$n_{p,c}(𝐫)=\frac{n(𝐫)T(𝐫)}{T_c(𝐫)},$$
(30)
where $`T`$ and $`T_c`$ are the temperatures of protons in the disc and the clouds, respectively. We take the temperature of the clouds to be the same everywhere - around $`T_c=10^5`$K (Kuncic et al, 1996), although the actual value will be absorbed into an overall constant - and find the temperature of the disc from the sound speed, $`c_s`$, as

$$T(𝐫)=\frac{m_pc_s^2(𝐫)}{k},$$
(31)
where $`k`$ is Boltzmann’s constant and $`m_p`$ is the mass of the proton. This gives us the cross section

$$A(𝐫)\left[\frac{kT_c}{m_pn(𝐫)c_s^2(𝐫)}\right]^{\frac{2}{3}}\left[\frac{1}{\rho (𝐫)c_s^2(𝐫)}\right]^{\frac{2}{3}}.$$
(32)
We take $`\rho (𝐫)`$ and $`c_s^2(𝐫)`$ from the solutions of the previous section. The expression should be normalised by the cross section of the clouds being around $`10^{10}`$ cm<sup>2</sup> at the innermost radius (Kuncic et al., 1996), however in practice we can absorb this constant into the constant $`\mathrm{\Gamma }`$, defined below.

The various assumptions above about the cloud structure (constant number of particles, constant temperature) result in clouds with size growing strongly with $`r`$, as $`r^{\frac{5}{3}}`$. Combined with the radial dependence of the flux emission (discussed in more detail below), this results in profiles dominated by the outer disc regions. We also want to consider more general radial dependencies, where profiles are dominated by the central regions. We therefore add the term $`r^g`$, to the cloud cross section, where generally $`g`$ is a small positive number. This has a somewhat analogous róle to the $`r^q`$ term in the standard thin disc models (Fabian et al. 1989), although has a physically different origin. The final expression is then

$$A(𝐫)\left[\frac{1}{\rho (𝐫)c_s^2(𝐫)}\right]^{\frac{2}{3}}r^g.$$
(33)
We take the number density of clouds to trace the number density of protons in the disc, hence

$$n_c(𝐫)=\mathrm{\Gamma }\rho (𝐫).$$
(34)
The constant of proportionality, $`\mathrm{\Gamma }`$, will determine the degree of shadowing by optically thick clouds. The probability of obscuration, $`P(𝐫)`$, with maximal value 1, is proportional to $`\mathrm{\Gamma }`$.

### 3.5 Emitted flux

We take the emitted flux to be proportional to the flux illuminating the clouds, $`f_{in}`$, the number density of the clouds, and the cross section of the clouds.

We take two contributions to the flux illuminating the clouds, a central illumination (e.g. due primarily to synchrotron radiation in the innermost regions of the flow for ADAFs Mahadevan 1997) an X-ray bath throughout the disc (due to Bremsstrahlung radiation). The central illumination falls of as $`\frac{1}{r^2}`$ and the Bremsstrahlung is proportional to the square of the number density of particles in the disc. So the emitted flux is given by

$`f_{out}(𝐫,𝐜(𝐫))n_c(𝐫)A(𝐫)f_{in}(𝐫)`$
$`\rho A(𝐫)^{\frac{2}{3}}\left[{\displaystyle \frac{A\mathrm{max}(\widehat{𝐫}\widehat{𝐜}(𝐫),0)}{r^2}}V_p(𝐫)+B\rho ^2\right].`$ (35)
Here $`A`$ and $`B`$ are constants determining the relative importance of central radiation and the X-ray bath. They will depend on the mass accretion rate of the system (Mahadevan 1997) and also the geometry of the clouds. The $`\mathrm{max}(\widehat{𝐫}\widehat{𝐜}(𝐫),0)`$ term contains the information that only the side of the cloud facing the centre of the disc is illuminated by the central flux. Hence, the corresponding reprocessed flux is only visible for clouds on the far side of the hole with respect to the observer, and the amount of flux seen by the observer decreases as the illuminated side of the clouds becomes more edge on to the observer. $`V_p(𝐫)`$ is the primary visibility function, analogous to the visibility function described earlier, giving the probability of photons emitted from the central regions reaching the cloud without being obscured by intermediate clouds. Fortunately, the photon paths from the central illumination are radial and so the integration of equation 15 becomes analytically tractable, using appropriate expressions for the cross sections and cloud number densities.

To give some visualisation of the processes involved, we will consider five cross sections of the disc at different heights above the equatorial plane, as shown in Fig.1. For each cross section, we compute a frequency map and a visibility map for both a centrally illuminated (synchrotron dominated) case and an X-ray bath (Bremsstrahlung dominated) case. The disc visualised is a self-similar ADAF disc with $`ϵ=0.1`$ and $`\alpha =0.35`$, viewed at 60 degrees. The results are shown in Fig.2, in which the observer is situated to the right.

Some thought reveals that the plots of Fig. 2 make sense physically, and the qualitative features deserve some further comments. The redshift map of the equatorial is different in two respects from the standard redshift maps of thin discs. Firstly, the regions of largest Doppler shift are rotated slightly with respect to the observer. This is because the radial motion is of similar magnitude to the rotational motion. Secondly, gravitational redshift is much more dominant. This is because the velocities are significantly less than the Keplerian velocities of a thin disc (see comments in §2.1). As we would expect, the redshifts become less pronounced as we move away from the equatorial plane.

The visibility in the centrally illuminated case is a superposition of the visibility in the X-ray bath case with the primary visibility. The primary visibility implies that the observer can only see clouds on the far side of the hole (hence there is more obscuration on the top part of the disc than the bottom part) and that the more edge-on clouds appear less illuminated.

The visibility in the X-ray bath case is determined by the size and number density of clouds as well as the the length of the of the trajectory through the clouds. This results in the central regions being the most obscured because the disc under consideration has $`g=1.5`$, which means the cloud number density times cross section is largest in the central regions. Also, the bottom regions are more obscured than the top regions because they are further from the observer.

## 4 Line profiles

Discs were considered $`(i)`$ at two or three different inclinations (15,30 and 60 degrees, or some subset thereof), $`(ii)`$ once with shadowing ($`\mathrm{\Gamma }`$ finite) and once without ($`\mathrm{\Gamma }=0`$), and $`(iii)`$ once with the profile dominated by the outer regions ($`g=0`$) and once with the profile dominated by inner regions ($`g=1.5`$). All discs have $`\alpha =0.35`$.

### 4.1 Self-similar ADAF profiles

We consider three ADAF discs. The first has $`ϵ=0.1`$, hence the density dependence is essentially spherically symmetric. We take this disc to be dominated by central illumination (i.e. $`A=10,B=1`$ in equation 3.5). The second disc also has $`ϵ=0.1`$ but is dominated by the X-ray bath (i.e. $`A=1,B=100`$). The third has $`ϵ=10`$, hence the density is concentrated within around 30 degrees of $`\theta `$ on each side of the equatorial plane. We take this disc to be dominated by central illumination.

A striking feature of the profiles from the two $`ϵ=0.1`$ ADAF discs considered, Figs. 3 and 4, is the lack of significant dependence on the inclination angle (compare to the profiles from flat thin discs in, e.g. Fabian et al. 1989, which show a strong inclination dependence). This is because the velocities given by equations 2-4 have magnitudes smaller than the Keplerian velocities of thin discs and so extreme blue Doppler shifts do not appear. There is more blue-shifted flux for the higher inclination discs, which is why the peak very near the rest frequency is lower in the higher inclination discs.

The discs with $`g=1.5`$ have profiles dominated by flux from the inner regions, which explains the greater presence of gravitationally redshifted flux. This is most notable in the X-ray bath case of Fig. 4 because here there is a contribution from Doppler redshifted flux which is not so visible in the centrally illuminated case. See Fig.2 for how the shadowings and redshifts interact. The fact that $`v_rv_\theta `$ in the discs of Figs. 3 & 4 is important because it rotates the most Doppler redshifted region into the less visible sector of the disc for the centrally illuminated case, Fig. 3.

The most notable feature of the profiles from the $`ϵ=10`$ ADAF disc, in Fig.5, is that cases with $`g=0`$ are of virtually identical form as the standard thin disc profiles, and show the same inclination dependence as the thin discs! This is not too surprising because the velocities of equations 2-4 give a magnitude very close the Keplerian magnitude and the angular motion is more important than the radial motion.

The profiles with $`g=1.5`$ are dominated by emission from the central regions, which results in the large amounts of red and blue shifted flux due to Doppler and gravitational effects.

### 4.2 Radial infall profiles

We consider a radially infalling disc dominated by central illumination. The profile is independent of inclination because of the spherical symmetry, so only one inclination angle is considered. The profiles are shown in Fig.6. To show explicitly the effect of shadowing, the non-shadowed profiles are plotted to scale with the shadowed profiles in the left column.

The profiles show a significant amount of blueshifted flux and little redshifted flux. This is the immediate consequence of the central illumination which means the blobs falling inwards towards the observer are visible, whilst those falling inwards away form the observer are illuminated on the side obscured from the observer. In the case where the profile is dominated by contributions from the central regions ($`g=1.5`$), we see a gravitational redshift and more blueshift from the fast inner flow.

### 4.3 Outflow profiles

We consider a disc with the bipolar outflow solutions of equations 7 to 11. The disc has $`ϵ=0.1`$ and is dominated by central illumination.

The profiles are very similar to those generated by the ADAF solutions without outflows. Here we show the 30 degree inclination profiles, but the profiles for different inclinations are very similar, as for the ADAF cases. The double peak for the outer region dominated case is seen in the top row of Fig. 7. The redshift dominates the flux in both cases. This is because of the gravitational redshift, and because the use of central illumination means that the observer primarily sees emission from the clouds moving away from his/her line of sight. This is opposite to the radial infall case.

## 5 Discussion and conclusion

We have considered a model in which cold, dense clouds are embedded in an otherwise optically thin, but geometrically thick accretion flow. Assuming these clouds reprocess primary X-ray emission from a central source and an X-ray bath, we have computed the expected fluorescent iron line profiles. This is a different scenario to the standard calculation of line profiles from thin, flat discs.

We considered various solutions for the dynamics of the disc: ADAFs, pure radial infall models, and bipolar outflows. We allow these solutions to specify the velocity, number density and size of the clouds at different points of the engine. In calculating the line profiles, a new effect compared to the thin disc models is the obscuration of flux by the clouds themselves, which are optically thick.

Many of the resulting line profiles are similar to those from thin disc models, suggesting that thin discs do not provide a unique model that is able to reproduce observations of iron line profiles. This suggests that cloud models warrant further study. Previous work has focused primarily on the X-ray continuum (Sivron and Tsuruta 1993; Bond and Matsuoka 1993; Nandra and George 1994; Collin-Souffrin et al. 1996) but there has been little modelling of the relativistic iron fluorescent line for the cloud model. Karas et al. 2000 considered non-relativistic profiles with large clouds at one radius.

A further incentive to consider thick disks with embedded clouds is that they may help resolve problems arising in thin disc modelling which relate to the inclination of discs and the reprocessed fraction. This includes the puzzling preference for low inclinations for Seyfert 2s (Turner et al. 1998), some incongruences between inferred inclinations from iron line modelling and those determined from other methods (Nishiura et al. 1998; Sulentic et al. 1998a; Wang et al. 1999), and the observed reprocessed fraction appearing to be too low and/or uncorrelated with line width (e.g. Lee et al. 2000; Chiang et al. 2000; Weaver 2000; Sulentic et al. 1998b). In this regard, we find that many of our embedded cloud, thick disc profiles (modelled with ADAF type flows) are relatively insensitive to the inclination angle because a significant pressure support means that the bulk velocities are smaller than the Keplerian values.

This latter feature also means that many of our profiles have sharp blue cut-offs because they are dominated by the gravitational redshift. Whether the sharp blue cutoff remains true in the Kerr case, where the disc extends closer to central hole and hence can acquire greater velocities should be investigated. Note that this will also depend fundamentally on the innermost radius at which the reprocessing clouds can survive. A thick geometry can also allow a wider range of reprocessed fractions, depending on the cloud distribution and amount of obscuration.

To summarise, the embedded cloud+ thick disc model allows a wider range of disc flow solutions to show reprocessed X-rays lines, and highlights how an alternative geometry can produce profiles similar to those from thin discs. The model can also lead to profiles with different properties to those of the thin discs, such as less sensitivity to inclination angle. As such, we suggest the cloud model deserves further investigation. This would include more consideration of the mechanisms of cloud formation and sizes, generalisation to the Kerr metric, and application to specific sources.

## Acknowledgements

Thanks to Clovis Peres for many fruitful discussions and early collaborations on this subject. SAH acknowledges the hospitality of the University of Rochester and the REU program. EB acknowledges partial support from DOE grant DE-FG02-00ER54600.


text:
# References

Version Date : 21th of July 1999

The broad Brillouin doublets and central peak of KTaO3

E. Farhia,b, A.K. Tagantsevc, B. Hehlena, R. Curratb, L.A. Boatnerd, and E. Courtensa

a Laboratoire des Verres, Université de Montpellier 2, F-34095 Montpellier, France

b Institut Laue Langevin, BP 156, 38042 Grenoble Cedex 9, France

c Laboratoire de Céramique, EPFL, CH-1015 Lausanne, Switzerland

d Solid State Division, Oak Ridge Nat. Lab., Oak Ridge, TN 37831-6056, USA

## Abstract

The incipient ferroelectric KTaO3 presents low-$`T`$ Brillouin spectra anomalies, *e.g.* a broad central peak (CP), and some additional Brillouin doublets (BD), whose origin is interpreted in terms of phonon-density fluctuation processes. A parameterisation from new extensive high-resolution neutron-scattering measurements is used to show that hydrodynamic second sound from high damping (compared to BD frequency) TA phonons may exist in the crystal. Furthermore, low damping thermal phonons may scatter light through two-phonon difference processes and appear on the Brillouin spectra either as a sharp or a broader BD, depending on the phonon damping and group velocity . The comparison between computed anisotropies and experimental measurements favours the second process.

keywords : Computer modelling, Perovskite compounds, Phonons, Spectroscopy.

Corresponding Author : E. Farhi

Address : Institut Laue-Langevin, BP 156, F-38042 Grenoble Cedex 9,France

Fax : 33 (0)4 76 48 39 06

E-mail Address : farhi@ill.fr

PACS numbers : 77.84 Dy, 63.20.-e, 62.20.Dc, 78.70.Nx, 66.90+r, 78.35+c

Physica B Preprint

Incipient ferroelectrics SrTiO3 or KTaO3 are perovskite type crystals that remain paraelectric down to the lowest temperature, due to zero-point quantum fluctuations . The dielectric constant saturates at a high value as the soft transverse optic (TO) zone-centre frequency stabilises to a low finite value. Since the case of SrTiO3 is complicated by the structural transition at $`T_a`$ 105 , we study KTaO3 which remains in the cubic phase down to the lowest temperatures.

Experimental evidence of anomalous behaviour in pure KTaO3 when cooling into the quantum paraelectric (QPE) regime are numerous. In Particular, the low energy Brillouin spectra present unusual features. Hehlen et al.  found a new sound-like Brillouin doublet (BD) both in SrTiO3 and KTaO3 in addition to the quasi-elastic broad central peak (CP) first reported by Lyons and Fleury . It was interpreted as the signature of second sound . Our purpose here is to investigate the origin of CP and BD in the frame of phonon-density fluctuation processes.

The anisotropy of anomalous Brillouin contributions was measured in KTaO3 Ultra High Purity crystals along the three high symmetry axes by means of Brillouin light-backscattering, with a momentum transfer $`|\stackrel{}{q}|=`$ 3.55 10-3 rlu. The observations of Hehlen et al.  are confirmed for the $`\stackrel{}{q}`$ // \<001\> symmetry axis. A sharp BD appears upon the broad CP  below $`T`$ 22 K. Its associated velocity is $`v_{BD}`$ 1100 m/s. In the \<110\> direction, the BD and CP contributions are also observed, as well as an additional broad BD (BBD) for temperatures between 30 and 100 K. While the CP and BBD contributions are also seen along the \<111\> direction, no sharp BD could be measured for temperatures as low as 5 K .

In order to look for an explanation of CP, BD, and BBD contributions in terms of phonon-density fluctuation processes, we parameterised the low-energy phonon dispersion surfaces $`\omega (\stackrel{}{q})`$ below 100 K within a simple phenomenological quasi-harmonic anisotropic model due to Vaks . The precise determination of the necessary parameters required additional new inelastic neutron-scattering experiments to be carried out. Taking into account the 4-dimensional nature of the instrumental resolution ellipsoid and the anisotropy of the dispersion surfaces, an effective energy uncertainty as low as 0.1 meV was achieved. The resulting dispersion curves and their parameterisation are presented in Fig. 1 along the principal symmetry axes, showing the well-known TO-TA coupling . The transverse phonon dispersion surfaces show deep *’valleys’* along the \<001\> directions and for $`\stackrel{}{q}`$ // \<110 \> with $`\stackrel{}{e}`$ // \<001\>. Some *’uphill’* regions are found elsewhere, specially for $`\stackrel{}{q}`$ // \<111\>, and $`\stackrel{}{q}`$ // \<110 \> with $`\stackrel{}{e}`$ // \<1$`\overline{1}`$0\>. This anisotropic character was confirmed by off-principal-axis measurements.

We now analyse the origin of the anomalous Brillouin contributions in terms of classical phonon-kinetics processes. Second sound (2S) is the collective mode of phonon-density fluctuations induced by multiphonon processes . It may propagate at a frequency $`\omega _{2S}`$ in a pure crystal when the ’window’ condition, $`\mathrm{\Gamma }_R<\omega _{2S}<\mathrm{\Gamma }_N`$ (1), is fulfilled. Here, $`\mathrm{\Gamma }_N`$ and $`\mathrm{\Gamma }_R`$ are the normal and resistive damping coefficients respectively. In KTaO3, Eq. (1) is satisfied since the proximity of the optic and acoustic modes enhances the relaxation processes. The computation of 2S velocity using Vaks model and the thermalized part of the phonon spectrum gives $`v_{2S}`$ = 1100 $`\pm `$ 200 m/s, in excellent agreement with measurements of BD velocity. However, it seems difficult to explain the directional dependence of BD in this framework.

When Eq. (1) is not satisfied, low-damping phonons might interact *via* two-phonon difference scattering (TPDS). Following Wehner and Klein , we evaluated the corresponding Brillouin spectra for different anisotropic dispersion surfaces. It was found that TPDS from low-$`\mathrm{\Gamma }_N`$ *valley*-TA with small group velocity $`v_g`$ appears as a sharp BD on low-$`T`$ spectra, whereas a BBD is obtained by TPDS at slightly higher $`T`$ from thermalized TA phonons around *uphill* regions, where $`\mathrm{\Gamma }_N`$ and $`v_g`$ are higher. Finally, a broad CP originates from high-$`\mathrm{\Gamma }_N`$ TO phonon TPDS. Thus, the intrinsically anisotropic TPDS spectral shapes are in agreement with Brillouin experiments.

As a conclusion, we reported new Brillouin spectroscopy and high-resolution inelastic neutron-scattering results in pure KTaO3. On the basis of those experimental data, a simple phenomenological parameterisation of the phonon dispersion surfaces was used to study two different scenarios for the appearance of additional Brillouin contributions : two-phonon difference scattering processes (TPDS) and second sound (2S). The comparison between computed anisotropies and experimental measurements seems to favour the TPDS explanation.