text:
# Ab initio and finite-temperature molecular dynamics studies of lattice resistance in tantalum

## I Introduction

The study of the plasticity of crystalline materials is a rich many-body problem involving physics on multiple length scales, with many remaining unexplained mysteries. The plasticity of bcc metals, for instance, is particularly challenging. Unlike their fcc and hcp counterparts, the bcc metals exhibit many active slip planes, have a strong temperature dependence in their plasticity and violate the simple empirical Schmid law schmid . Moreover, theoretical calculations of the most basic question in plasticity, the stress needed to induce yield at low temperature in a pure sample, differ from experimental extrapolations by over a factor of two duesbery3 . The purpose of this work is to provide needed insight into this discrepancy.

Ultimately, it is the physics of the $`111`$ screw dislocation defect which controls the low-temperature plasticity of bcc materials vitek ; duesbery3 . The Peierls stress, the yield stress at which these dislocations first begin to move spontaneously, is difficult to compare directly with experiment. Whereas most computational work on the Peierls stress measures the stress to move an isolated, infinitely straight dislocation at zero temperature vitek ; duesbery1 ; segall ; wang ; rao1 ; woodward1 ; woodward2 ; yang ; xu ; moriarty ; mgpt ; takeuchi ; ito , experiments measure the Peierls stress at a finite temperature in systems with many interacting, curved dislocations and in media with defects and surfaces. As an example of the present challenges which the literature faces, using model generalized pseudopotential theory (MGPT) mgpt , Yang et. al.yang  predict for the T=0 Peierls stress a value 2.5 times greater than experimental extrapolations tang . Because such potentials are not based upon first principles, it is impossible to determine a priori whether this discrepancy is due to the inter-atomic potential, the environmental complexities discussed above, or to a flaw in our understanding of the relation between the Peierls stress and the experiments.

As it is a daunting experimental task sigle  to observe properties of a single dislocation locked deep in the heart of a material, accurate theoretical calculations of such systems is essential. Nearly all theoretical calculations to date, concerning such dislocations, have relied upon empirical potentials vitek ; duesbery1 ; segall ; wang ; shenoy ; moriarty ; xu ; mgpt ; ito ; yang ; rao1 ; bulatov ; takeuchi . Given the empirical nature of such calculations, the complex directional bonding properties of bcc materials, and the lack of direct comparison with experiments for validation, first principles ab initio calculations of dislocations in such systems are clearly needed. Ismail-Beigi and Ariassohrab\_disloc  were the first to show that density functional theory calculations were crucial in understanding the fundamental properties of the $`111`$ screw dislocation core structure in bcc molybdenum and tantalum. Until that work, most computational studies based on empirical potentials vitek ; duesbery1 ; rao1 ; mgpt ; wang ; xu ; moriarty ; segall , supported the idea that the dislocation core breaks symmetry, with two energetically equivalent ground state structures which spread outward along two different equivalent sets of three $`\{110\}`$ planesvitek , similar to the concept originally proposed by Hirsch hirsch ; mitchell . Until the availability of the ab initio calculations, the prevailing view of the violation of the Schmid law in the bcc metals was based upon this structure vitek . Ismail-Beigi and Arias sohrab\_disloc , in contrast, showed that for both molybdenum and tantalum the ground-state structure within density functional theory was a non-degenerate symmetric core, strongly supporting the work of Suzuki and Takeuchi suzuki ; takeuchi ; kimura  which first suggested that it is the Peierls potential itself that controls the lattice resistance and not the details of the core structure. To help resolve the discrepancy between theoretical and experimental Peierls stresses, the work below provides a reliable ab initio prediction of the Peierls stress in bcc tantalum which is free of the unrealistic electronic boundary conditions employed in the only other ab initio prediction of the Peierls stresswoodward1 ; woodward2 .

Here, we show that the Peierls stress, calculated within density function theory, is over a factor of five larger than expected from extrapolation of experimental results tang . This supports the view that the discrepancy between the experimental and computational predictions are largely due to the aforementioned environmental complexities, to a flaw in relating the experimental data to the Peierls stress, or to a combination of both.

To further explore possible physical effects leading to this discrepancy, we study the extrapolation of experimental data to determine the zero-temperature Peierls stress. Such extrapolations generally employ fits from mesoscopic or thermodynamics/kinetic models tang ; kocks ; cuitino ; stainer . However, it has not been established that such models can accurately describe the lowest temperature regime correctly, placing doubt on the quality of these extrapolations. To address this issue, the work below also provides the first temperature- and orientation- dependent study of the Peierls stress in a bcc metal. We moreover show that extrapolation of our finite temperature results using a current fitting model leads to an underestimation of the zero-temperature Peierls stress. This underscores the difficulty in extrapolating the experimental data accurately but does not fully account for the observed discrepancy.

In Section II, this manuscript reviews the various techniques in use for calculation of the Peierls stress in the context of efficacy for application to ab initio calculations. Section III gives the first calculation of the temperature and orientation dependent Peierls stress in a bcc material. Section IV describes our technique for obtaining Peierls stresses within small unit cells with periodic boundary conditions. Finally, Section V presents our ab initio prediction for the Peierls stress and compares and contrasts it to currently available experimental and computational values.

## II Boundary Conditions

The fundamental distinction among theoretical approaches to calculation of the Peierls stress is the choice of boundary condition. The literature describes three types of boundary conditions: cylindrical boundary conditions vitek ; segall ; shenoy ; mgpt ; xu ; takeuchi ; ito ; duesbery1 , Greens function (or “flexible”) boundary conditions sinclair ; rao2 ; yang ; rao1 ; vitek ; woodward1 ; woodward2 , and periodic boundary conditions segall ; wang ; bulatov . We now briefly review each with emphasis on the unique challenges of ab initio electronic structure calculations.

### II.1 Cylindrical Boundary Conditions

In the practice of cylindrical boundary conditions, anisotropic elasticity theory hirth ; stroh ; head  is used to generate a dislocation in the center of a cylinder. The cylinder is then separated into inner and outer regions. The atoms in the outer region are held fixed to the solution of anisotropic elasticity theory while the atoms in the inner region relax under the inter-atomic forces. To calculate the Peierls stress, a stress is applied to the system until the dislocation moves.

This approach suffers numerous draw backs when applied to density functional theory. To avoid surface effects and to properly account for the non-linear nature of the dislocation, such cylinders generally have to be quite large. First, even the outer cylinder is of finite size and therefore the outer region must be sufficiently large enough so that forces generated by it onto the inner region are equivalent to those generated from an infinite continuum. The inner region also must be sufficiently large to mitigate two effects. The inner region must be large enough so that linear elasticity theory represents well the forces which it imposes on the outer region. The inner region also must be large enough so that motion of the dislocation is not adversely affected by the fixed outer region, which is a concern because the fixed outer region represents the displacement field when a dislocation is at its center and therefore generates a extraneous force which tends to prevent motion of the dislocation shenoy . When using simple, inter-atomic potentials, the use of large cylinders mitigates all of these effects. However, this approach is not viable for density functional calculations with their extreme computational demands.

This approach, moreover, is particularly ill-suited for electronic structure calculations because the artificial surface at the outside of the outer region, being far different from the bulk, give rise to strong scattering of the electrons far different than would an infinite continuum. This is particularly problematic for metals, because the single-particle density matrix, which quantifies the effects of this scattering on the inter-atomic forces, decays only algebraically in metals sohrab\_densmat . The following subsection demonstrates that the boundary regions should be quite large in order for these surface effects not to result in large fictitious forces in the active region of the calculation.

### II.2 Greens Function Boundary Conditions

The use of Greens function, or flexible, boundary conditions sinclair ; rao2 ; vitek , is an effective way to reduce the size of the simulation cell. This approach also employs a cylindrical geometry. However, rather than the “inner” and “outer” atomistic regions of the cylindrical boundary approach, the Greens function approach employs three inter-atomic regions: an inner “core” region containing the center of the dislocation, an intermediate “buffer” region, and an outermost “continuum-response” region. With proper implementation, the outer and inner regions couple only indirectly through the response of the buffer region.

In this method, all three regions respond to the presence of a dislocation; however the response of each region is treated differently through a number of steps. Initially, all regions are displaced by the solution to anisotropic elasticity theory. Each iteration then begins by relaxing the atoms in the core region according to the forces which they experience, as computed from either an inter-atomic potential or an ab initio method. The forces generated from the mismatch between the outer and inner regions, which the cylindrical approach above ignores, are then relieved by moving the atoms of all three regions according to the elastic Greens function solution, leaving only the nonlinear effects from the core region unaccounted. The next iteration then begins by relaxing these forces as described above. Iterations proceed until until the forces in the core and buffer region are negligible.

What distinguishes this approach from simple cylindrical boundary conditions is that the continuum region, via the Greens function response, is allowed to respond to the motion of the dislocation and to the elastic response generated by the core region as the dislocation moves. So long as the continuum response region (a) accurately represents the structure induced by the presence of the dislocation and (b) is sufficiently wide to properly reproduce the forces on the atoms in the buffer and inner regions, this approach accurately describe basic properties of a dislocation.

In order for the first assumption (a) above to hold, the inner core region must be sufficiently large to contain all atoms with displacements outside of the linear regime and the buffer region must be sufficiently wide so that displaced atoms in the core have no effect on the forces experienced in the continuum-response region. The second assumption (b) requires that the continuum-response region to be sufficiently large so that its termination has no effect on the forces on the atoms in the buffer or inner region. The radius of the calculation must therefore exceed the sum of the non-linear core radius plus twice the range over which motion of atoms creates forces within the lattice. As the latter range can be quite large for electronic structure calculations in metals, the application of this approach to electronic structure calculations can be problematic.

The Greens function approach has predicted successfully dislocation properties when applied to time consuming empirical potentials rao1 ; yang  which have a limited interaction range. The approach also has been applied to density functional calculations of the Peierls stress for molybdenum and tantalumwoodward1 ; woodward2 , where its application is more questionable due to the above interactions. In these latter works, the artificial boundary on the outside of the continuum-response region have been treated in one of two ways woodward1 ; woodward2 , either by keeping the surface free in vacuum or by embedding in periodic boundary conditions with the vacuum filled with material which must contain severe domain boundaries due to the incompatibility of a net Burgers vector with periodic boundary conditions.

To gauge the effects which this artificial boundary may have and how far these effects penetrate from the continuum-response region into the buffer region, we perform a test calculation within the density functional theory pseudopotential approach payne  of the magnitude of the forces generated onto the system due to the presence of a domain boundary similar to those in the works cited above woodward1 ; woodward2 .

For this calculation, we employ the same computational procedure as for our production calculations in Section V. Here, however, as this is a test, we employ only a single $`k`$-point to sample the Brillouin zone ($`\mathrm{\Gamma }`$). We begin with an orthorhombic cell of 24 atoms of tantalum in a bulk arrangement with supercell lattice vectors $`\stackrel{}{r}_1=a[1\overline{1}0]`$, $`\stackrel{}{r}_2=4a[11\overline{2}]`$ and $`\stackrel{}{r}_3=a/2[111]`$. We choose this cell because its length along $`\stackrel{}{r}_2`$ is the same as the smallest simulation cell used in References woodward1,  and woodward2, . We then generate a domain boundary at the edge of cell along the $`(11\overline{2})`$ plane by changing the lattice vector $`\stackrel{}{r}_2`$ to $`\stackrel{}{r}_2=4a[11\overline{2}]+\alpha \stackrel{}{r}_3`$ and holding the atoms in the unit cell fixed in their bulk locations. $`\alpha `$ is chosen such that the shift is small and the nearest neighbor distance is always within $`95\%`$ of the bulk, representing even less of a disturbance that in Reference woodward1, , where atom were within $`90\%`$ of the bulk nearest neighbor distance. To estimate the effect of the scattering of electrons at the domain boundary on the inter-atomic forces, we hold the atoms fixed and compute the ab initio forces acting upon them.

Figure 1 shows the forces along the $`[111]`$ direction as a function of distance from the center of each domain. Note that relatively large forces develop deep within the cell. This data indicates that the continuum-response region should be quite large ($`510`$ Å) in order for the response of the electrons not to adversely effect the forces in the buffer region. Note also that the the buffer region should be of similar width to prevent forces from the non-linear displacements in the core from penetrating into the linear continuum-response region. Such large continuum-response and buffer regions can make the calculation infeasible with current computational techniques.

In fact, the only density functional calculations of the Peierls stress in this system to date employ the Greens function method but with a distance from the buffer region to the domain boundary of only $`3.7`$ Å. It thus is unclear whether the continuum region in these calculations is sufficiently large to lead to reliable results and clearly further calculations are needed to support those results. Below, we provide just such calculations using the method of periodic boundary conditions, which perturb the electronic system far less than the introduction of domain boundaries.

### II.3 Periodic Boundary Conditions

The final common choice for boundary conditions is to repeat the dislocation core periodically throughout space, so that the dislocation is no longer isolated, but embedded in bulk material containing an array of dislocations. Consistency with periodic boundary conditions then demands that the unit cell contain an even number of dislocations with Burgers vectors of alternating sign, arranged typically in dipolar marklund ; ariasjoannop ; lehto  or quadrupolar bigger ; lehto ; sohrab\_disloc  arrays. For static properties of the dislocation core, such cells give reliable results as the elastic fields of the surrounding dislocations effectively cancel at the location of each core.

Periodic boundary conditions can also been used to calculate the Peierls stress segall ; wang ; bulatov . Care must be taken, however, because the symmetry of the dipole/quadrupole array breaks as the dislocations move and the dislocation-dislocation interaction then ceases to be negligible. The use of large unit cells can control this effectbulatov ; however, such a direct, brute-force approach is not practical for computationally demanding ab initio calculations. To make such calculations feasible, they must occur in small periodic cells, thereby demanding proper accounting of the dislocation-dislocation interactions.

We have shown in another work segall  that, under certain conditions, such interactions can be accounted accurately with minimal extra computational effort, so that accurate values of the Peierls stress can be obtained from density functional calculations in periodic cells. As the residual errors with this approach are associated with the boundary conditions, the magnitude of such an error can be tested by other computational methods, such as empirical potentials. Our previous work shows that that this residual error is relatively small segall , a fact which we confirm explicitly below.

Because the deviations from the bulk arrangement at periodic boundaries are relatively mild, such calculations are ideal for mitigating electronic boundary effects. Given the simplicity of working with these boundary conditions and the possibility of the extraction of accurate values for the Peierls stress from small unit cells, we choose to work with periodic boundary conditions throughout this work. Section IV outlines our procedure for calculating the Peierls stress while working with periodic boundary conditions and describes the sources and the magnitude of the residual errors. (See Reference segall,  for a full discussion of these issues.)

## III Dependence of the Peierls stress on orientation and temperature

To illustrate the complexities of relating computational predictions to experimental findings, we now explore the dependence of the Peierls stress in bcc tantalum on orientation and temperature. The strong dependencies which we shall find underscore the unique properties of dislocations in bcc metals. To clarify, as some authors use slightly different definitions for the Peierls stress, here we consider the Peierls stress as the value of the stress on the maximum resolved shear stress plane (maximum value of the shear stress along the -direction) when the dislocation first moves to a different equilibrium position.

Despite recent advances in ab initio quantum mechanical methods, such methods are still too computationally intensive to study such properties as the temperature dependence of the Peierls stress. Therefore, for these calculations, we employ a molecular dynamics (MD) framework carried out using a first-principles-based, many body force field (FF) for tantalum, which we denote qEAM, which we have developed to allow accurate and computationally efficient evaluation of atomic interactions strachan ; wang .

As described above, we carry out these calculations within periodic boundary conditions. The super-cell consists of a quadrupolar arrangement bigger ; lehto  of dislocations containing 5670 atoms with lattice parameters $`a_x=70.59`$ Å, $`a_y=73.39`$ Å and $`a_z=20.11`$ Å, where the $`x`$-, $`y`$-, and $`z`$\- axes of our coordinate system are along $`[1\overline{1}0]`$, $`[11\overline{2}]`$, and $`[111]`$ directions, respectively. As we have shown in another work segall , such a cell gives very accurate values for the Peierls stress.

### III.1 Orientation dependence of the zero-temperature Peierls stress of the $`111`$ screw dislocation

To calculate the zero temperature Peierls stress, we start with a fully relaxed quadrupole dislocation configuration at zero stress and increase the stress in steps of 50 MPa until the dislocations move. Once the dislocations move, we restart the calculation from the structure equilibrated just prior to the motion and increase the stress in smaller steps (5 MPa) in order to more narrowly define the critical stress. At each incremental target stress, we relax the atoms and stresses in the cell by running two very low temperature (T=0.001 K) MD simulations. The first run is for 15 ps at constant stress and temperature (N$`\sigma `$T ensemble) using a Rahman-Parrinello barostat rahman  and a Hoover hoover  thermostat, and the second run is for 50 ps at constant volume and temperature (NVT ensemble). We find this approach to be quite stable for relaxing the cell and the atoms of the system.

The $`111`$ screw dislocation has three equivalent {112} and three equivalent {110} potential slip planes, with such planes occurring at 30 intervals. To study the orientation dependence of the Peierls stress, we apply three types of pure shear stress to the system: a $`\sigma _{xz}`$ stress, a positive $`\sigma _{yz}`$ stress and a negative $`\sigma _{yz}`$ stress. (Note that with coordinate axes as defined above, the $`z`$-axis lies along the dislocation line.) These stresses leads to forces on the dislocation in the $`112`$, $`110`$-twinning and $`110`$-antitwinning directions, respectivelyhirth . Along these directions, we find Peierls stresses of $`\tau _{112}`$ = 655 MPa, $`\tau _{twin}`$ = 575 MPa, and $`\tau _{antitwin}`$ = 1075 MPa, respectively. Table 1 shows that our essentially zero-temperature results are in good agreement with those of Yang and collaborators yang , who employed model generalized pseudopotential theory (MGPT)mgpt , a different inter-atomic potential. The result of such a strong dependence of the Peierls stress on orientation is consistent with the experimentally observed breakdown of the Schmid law in bcc metals.

To make quantitative comparison with experiments, which are carried out at nonzero temperature, we compare our results to those of Tang et. al. tang , who fit experimental data wasserbach  to a mesoscopic model and then extrapolate to extract the zero-temperature Peierls stress. Their predicted value of 248 MPa for the $`112`$ Peierls stress is over a factor of two lower than our result. This type of discrepancy, where the theoretical Peierls stress overestimates the zero-temperature extrapolation of the experimental data by a factor of two to three, is quite generally observedduesbery3 . This discrepancy may be due either to inaccuracies in the theoretical calculations or, perhaps, to a flaw in the comparison between the zero-temperature extrapolation of the experimental data and theoretical predictions.

### III.2 Temperature dependence of the Peierls stress of the $`111`$ screw dislocation

To explore potential difficulties with the zero-temperature extrapolation, we now present what to our knowledge is the first temperature-dependent study of Peierls stress using a realistic potential for a bcc metal. For these calculations, we continue to employ the qEAM FF and begin with the zero-temperature, equilibrated structures. We then apply various constant shear stresses (lower than the T=0 Peierls stress) to the system while slowly increasing the temperature (in steps of 10 K) until the dislocations move. Similarly to the T = 0.001K case, for each temperature we first run for 10 ps in the N$`\sigma `$T ensemble and then for 25 ps in NVT ensemble.

Because the Peierls stress can depend on the rate at which the strain is applied, to place our results in context, we first estimate the strain rate in our computations. The strain rate is approximately $`\dot{\gamma }=\rho v_db`$, where $`v_d`$ is the dislocation velocity, $`\rho `$ is the dislocation density and $`b`$ is the Burgers vector. Using a dislocation density typical of the experimentstang  ($`\rho =10^{11}1/m^2`$) and estimating the dislocation velocity as the ratio between the distance traveled in one jump ($`1/3a112=2.717`$ Å) and the simulation time (35 ps), we obtain an effective strain rate of $`10^21/s`$, which is large compared to the strain rates ($`4\times 10^5`$) in the experiments used for the zero-temperature extrapolationstang ; wasserbach .

Figure 2 summarizes our results for the temperature dependence of the Peierls stress as a function of temperature for the three directions ($`112`$, twinning and anti-twinning). As expected, the Peierls stress obtained from our MD simulations decreases rapidly with increasing temperature, particularly for very low temperatures. It is important to mention that, although our simulations are three dimensional, the dislocations move as straight lines without the formation of double kinks because our simulation cell is only seven Burger’s vectors long along the dislocation lines. Such double kinks are quite important at finite temperatures as they tend lower the lattice resistance at non-zero temperatures. Our results are approximately a factor of $`24`$ larger than the fit of Tang et. al. tang  to the experimental data of Wasserbäch wasserbach . We feel that this is reasonable, considering the facts that our simulation cells do not allow for double kink formation and that, as discussed above, our strain rates are much higher than those in the experimentstang ; wasserbach .

Experimental extrapolations of the Peierls stress to zero temperature generally come from mesoscopic or kinetic/thermodynamic model tang ; kocks ; cuitino ; stainer  fits to experimental data and extrapolated to zero temperature. To explore the effects of this procedure, we fit our atomistic data to such a model, perform the extrapolation and then compare with our direct zero-temperature results.

For the fit we use an analytical expression for the dependence of the Peierls stress with temperature at constant strain rate stainer ; cuitino  based on a mechanism involving double kink nucleation and propagation. This model gives for the temperature dependent Peierls stress

$$\tau _P=\frac{\tau _0}{\beta E^{\mathrm{kink}}}\mathrm{asinh}\left(\frac{\dot{\gamma }}{\dot{\gamma }_0^{\mathrm{kink}}}\mathrm{e}^{\beta E^{\mathrm{kink}}}\right),$$
(1)
where $`\beta `$ is $`1/k_BT`$, with $`k_B`$ being Boltzmann’s constant and $`T`$ the absolute temperature, $`E^{\mathrm{kink}}`$ is the kink energy hirth , $`\tau _0`$ is the effective Peierls stress and $`\dot{\gamma }_0^{\mathrm{kink}}`$ is the reference strain rate. Here, the effective Peierls stress is

$$\tau _0=\frac{E^{\mathrm{kink}}}{bL^{\mathrm{kink}}l_P}$$
(2)
and the reference strain rate is

$$\dot{\gamma }_0^{\mathrm{kink}}=2b\rho l_P\nu _D,$$
(3)
where $`b`$ is the Burgers vector, $`L^{\mathrm{kink}}`$ is the kink length, $`\rho `$ is the dislocation density and $`\nu _D`$ is the attempt frequency which may be identified with the Debye frequency to a first approximation stainer ; cuitino , and $`l_P`$ is the distance between two consecutive Peierls valleys. Physically, $`E^{\mathrm{kink}}`$ is the minimum energy to form a double kink, $`L^{\mathrm{kink}}`$ is the minimum length for this double kink and $`\tau _o`$ is the stress, whose work to move a dislocation a distance $`l_P`$ is equal to $`E^{\mathrm{kink}}`$.

Figure 2 shows the fit of Equation (1) to our atomistic data. To mimic how zero-temperature lattice resistances are generally extracted, we adjust the three unknown parameters ($`\tau _0,E^{\mathrm{kink}}\text{ and }\dot{\gamma }_0^{\mathrm{kink}}`$) to fit our higher temperature data only. Intriguingly, extrapolation of our higher temperature results to zero-temperature leads to an underestimation of the Peierls stress of between $`10\%`$ and $`20\%`$. We would also expect that a fit to data from a cell sufficiently large to allow for double kink formation (which is active in the experiments at non-zero temperatures) would lead to an even a larger underestimation of the zero-temperature stress. These results therefore suggest that the general discrepancy between extrapolated experimental values and the calculated values for the zero-temperature Peierls stress may be the result of failure of nonzero-temperature models to describe properly the low-temperature regime. This illustrates one possible difficulty in relating the experiments to computational predictions and underscores the need for first principles studies of such a phenomena.

## IV Accurate Peierls Stress Calculations in small periodic cells

Having underscored the need for first principles electronic structure studies and already determined the most effective boundary conditions for such studies as periodic, we now focus on determination of the cell of minimal size appropriate to calculation of the zero-temperature Peierls stress for a $`111`$ screw dislocation in a bcc metal when the maximum resolved shear stress is along a {110} plane.

To minimize image effects, we employ periodic boundary conditions with a quadrupolar unit cell. Figure 3 shows a differential displacement map vitek  of such a cell of size 42Å$`\times `$ 41Å in the plane perpendicular to the Burgers vector. In such a map, the dots indicate columns of atoms along the $`[111]`$, and the vectors between the columns of atoms indicate the relative shift along the Burgers vector due to the presence of a dislocation between each pair of columns, with the vectors scaled so that a vector of full length between the columns corresponds to $`1/3`$ a Burgers vector. In this ground state structure of the dislocation, triads of full length vectors surrounds the center of each dislocation, corresponding to a new displacement by a full Burgers vector upon completion of a closed loop about the each center. In practice, because of symmetry, the quadrupole cell may be reduced in half containing two dislocations, when lattice vectors are properly chosen.

As Section II notes, for calculations of static properties such as the ground-state dislocation core structure, the strain fields from the surrounding dislocations in a quadrupolar array essentially cancel at each dislocation core. However in a dynamical problem such as the calculation of the Peierls stress, the dislocations begin to interact with the stress fields of the others as they begin to move. Our previous worksegall  shows that, for the particular geometry considered here, accurate values for the Peierls stress can be extracted from quite small unit cells provided the proper procedure is followed. We now outline that procedure while reviewing the relevant background.

### IV.1 Calculation of Peierls stress within periodic boundary conditions

To calculate the Peierls stress in a small periodic cell, we begin with lattice vectors appropriate to bulk material in the absence of dislocations and then, while relaxing the internal coordinates of the cell, apply increasing pure $`ϵ_{xz}`$ strains (same Cartesian coordinates from Section II.3) until the dislocations move. Such a strain drives the dislocation along the $`[11\overline{2}]`$ directionhirth . Before the dislocation moves, the strain energy of the cell increases quadratically,

$$E=\frac{1}{2}C^{}ϵ_{xz}^2,$$
(4)
where $`C^{}`$ is an elastic constant associated with the quadrupole unit cell which can be extracted simply from the energies of the cell as the strain increases. The stress associated with this strain is

$$\sigma _{xz}=C^{}ϵ_{xz},$$
(5)
so that at the strain at which the dislocation moves, Equation (5) gives the Peierls stress.

The great benefit of the above procedure is that it requires a minimal search through phase space in order to calculate the Peierls stress and accounts accurately for the effects of the dislocation-dislocation interactions. The $`C^{}`$ elastic constant, which comes as a direct byproduct of the procedure and requires no further calculations, suffices to account accurately for the leading order effects of the inter-dislocation interactions. This correction captures most of the effects from working with small unit cells and can differ from that of an equivalent bulk cell by a value of two or more.

The above procedure involves several approximations requiring justification. First, we apply strain relative to the lattice vectors of the bulk cell in the absence of the dislocation, rather than those of the quadrupole array. Previously, we have shown through explicit calculation on model inter-atomic potentials that working with the relaxed lattice vectors of the quadrupolar array does not improve the value calculated for the Peierls stress and therefore is not neededsegall . The reason for this is that although working with the bulk lattice vectors generates artificial stresses, these are primarily diagonal ($`\sigma _{ii}`$), because the greatest effect of the presence of the dislocations is to dilate the systemsegall . In the present geometry, such diagonal stresses result only in a constant shift in the energy over the range of applied strain and do not generate driving (Peach-Kohlerhirth ) forces on the dislocations.

Some slight care must be taken with the above argument. Duesbery and others duesbery1 ; ito  have pointed out that stresses which do not result in driving forces on a dislocation still may affect the overall value of the Peierls stress needed to drive the dislocation, because such stresses may modify the dislocation core structure duesbery2 , an effect not accounted in linear elasticity theory. This effect of non-driving stresses is, in fact, one of the common violations of the Schmid law which bcc metals exhibit. Diagonal stresses, however, do not have a large effect duesbery1  on the value for the Peierls stress, as they tend to only compress or expand the core. Fortunately, because the bulk lattice vectors should be relatively close to the that of the dislocation cell, we expect all of these effects to be quite small, as we have found previouslysegall  and again verify in the test calculations below.

The second simplification in our procedure is that rather than applying a stain which imposes a pure $`\sigma _{xz}`$ stress, we apply a pure $`ϵ_{xz}`$ stain, which also generates a residual $`\sigma _{xy}`$ stresshirth . To generate a pure stress of the form $`\sigma _{xz}`$, one would have to apply an additional $`ϵ_{xy}`$ stain of a magnitude determined by yet another elastic constant of the quadrupolar array. Because the calculation of this constant would increase significantly the number of calculations required and because the residual $`\sigma _{xy}`$ stresshirth  acts on the plane perpendicular to the dislocation, and thus does not create a driving force on the dislocations, we simply apply the pure $`ϵ_{xz}`$ strain. As with the diagonal stress components, although the residual in-plane $`\sigma _{xy}`$ stress does not drive the dislocations, it can affect the Peierls stress by modifying the core structure. Unlike the diagonal stress components, the in-plane stress does significantly affect the Peierls stress in bcc metalsduesbery1 ; ito . This effect, however, will be small so long as the residual $`\sigma _{xy}`$ stress is small compared to the driving $`\sigma _{xz}`$ stress. The ratio $`\sigma _{xy}/\sigma _{xz}`$, is equal to $`C^{\prime \prime }/C^{}`$ where $`C^{}`$ is the elastic constant appearing in Equation (4) and $`C^{\prime \prime }`$ is another combination of elastic constants. In pure bulk cubic materials these constants have the formhirth

$`C^{}`$ $`=`$ $`{\displaystyle \frac{1}{3}}(C_{11}+C_{44}C_{12})`$
$`C^{\prime \prime }`$ $`=`$ $`{\displaystyle \frac{\sqrt{2}}{6}}(2C_{44}+C_{12}C_{11}),`$
where the $`C_{ij}`$ are the standard elastic constants for cubic materials.

As evidence of the correlation between the ratio of these constants and the errors in the Peierls stress, we note that in our previous studysegall  with same empirical potential as in Section III, the ratio $`C^{\prime \prime }/C^{}`$ varied from approximately $`1/5`$ in the smallest cell studied to less than $`1/10`$ for all other cells while the error in extracting the Peierls stress went from 18% to less than 2%, respectively. For the potential used in that study, the value of $`C^{\prime \prime }/C^{}`$ computed from bulk elastic constants is $`1/90`$. As the bulk value of $`C^{\prime \prime }/C^{}`$ within density functional theory is less than $`1/100`$sohrab\_disloc , we expect the errors in our ab initio value of the Peierls stress to be even somewhat smaller. Lending further support to this view is the result of Duesbery and Vitek duesbery2  demonstrating that the effect of $`\sigma _{xy}`$ stresses on core structure is much less for the non-degenerate core structure, which we have in our density functional calculations, than for the degenerate core structures, which we had in our inter-atomic potential calculations segall .

### IV.2 Demonstration

To demonstrate the efficacy of the above procedure, we now proceed to extract the Peierls stress in the $`112`$-direction for vanadium and tantalum from calculations in small periodic cells, when the maximum resolved shear stress is along a {110}-plane. In all calculations we calculate the Peierls stress for an infinite straight dislocation as our periodic cells has lattice vector $`\stackrel{}{a}_3`$ = $`a[111]/2`$ along the dislocation line. To allow comparison with the Peierls stress of isolated dislocations, we employ empirical potentials for this demonstration. For vanadium we use the Finnis-Sinclair finnis  potential with modifications made Ackland and Thetford ackland , and for tantalum we use the same potential as in Section III but with a slight adjustment of parameters to produce a non-degenerate core structure. The bulk ratios for $`C^{\prime \prime }/C^{}`$ for vanadium and tantalum within these models are $`1/10`$ and $`1/6.5`$, respectively, much larger than the density functional theory value.

To determine the reference value for the Peierls stress, we employed cylindrical boundary conditions with large amounts of material, increasing the radius of the cylinders until the boundary forces were small shenoy  and the Peierls stress approached an asymptotic value. To extract the Peierls stress from within periodic boundary conditions, we follow precisely the procedure which Section IV.1 outlines.

Figure 4 shows the resulting energy versus stain curve for vanadium within a periodic cell of size 42Å$`\times `$ 20Å. At a strain $`0.054`$, the curve exhibits a discontinuity signaling the critical strain for moving the dislocation. The curvature of the fit determines the elastic constant $`C^{}`$ through Equation (4). Finally, combining this value of $`C^{}`$ with the observed critical strain Equation (5) yields the Peierls stress. We repeated this procedure for tantalum as well.

Table 2 summarizes our results for both vanadium and tantalum. The table shows that the errors are relatively small, much smaller than the general discrepancy between empirical potentials and the experimental extrapolations. It is also noteworthy that the potentials employed in this demonstration exhibit $`C^{\prime \prime }/C^{}`$ ratios over an order of magnitude larger than density functional theory. From these and previous results segall , we conservatively estimate that the error in the Peierls stress in the density functional calculations below should be no greater than $`20\%`$.

## V Density functional results and discussion

### V.1 Computational Details

All of our first principles electronic structure calculations employ the plane-wave density-functional pseudopotential approach payne  within the local density approximation ceperley ; perdew . We employ pseudopotential of the Kleinman-Bylander form kleinman  with $`s`$, $`p`$ and $`d`$ non-local channels which has been used successfully in previous workswoodward\_ps ; sohrab\_disloc  and a plane wave basis with cutoff of 40 Rydberg. As justified above, we employ a super-cell containing a quadrupolar array of dislocations of size $`\stackrel{}{r}_1=5a[1,\overline{1},0],\stackrel{}{r}_2=(3/2)a[1,1,\overline{2}],\text{and}\stackrel{}{r}_3=a[1,1,1]/2`$, where $`a=3.25`$ Å is the lattice constant of the cubic unit cell. The lattice vectors of this cell are $`\stackrel{}{a}_1=\stackrel{}{r}_1/2\stackrel{}{r}_2+\stackrel{}{r}_3/2,\stackrel{}{a}_2=\stackrel{}{r}_1/2+\stackrel{}{r}_2+\stackrel{}{r}_3/2\text{ and }\stackrel{}{a}_3=\stackrel{}{r}_3`$. To carry out the integrations over the Brillouin zone we use a non-zero electronic temperature of $`k_BT=0.1`$ eV to facilitate integration over the Fermi surface and sample the zone at sixteen special $`k`$-points kpts . These choices give energy differences reliably to within 0.1 eV/atom. Finally, to determine the electronic structure, we minimize using the analytically continued functional approach arias , expressed within the DFT++ dftpp  formalism.

### V.2 Energy landscapes

To demonstrate the discrepancies that occurs between the predictions of empirical potentials and first principles electronic structure studies, we compare the energy landscape for a dislocation moving along a reaction coordinate from the easy core configuration to the hard core configurationvitek ; xu  for both molybdenum and tantalum as determined within MGPTmgpt  and as calculated ab initio. Within MGPT, we have carried out the calculation for molybdenum ourselves and we used the hard-easy core energy difference reported in the literaturemoriarty  for tantalum. Within density functional theory, we have calculated the energy at a number of points along the reaction pathway for tantalum, and for molybdenum we report the difference between the easy core and hard core configurations as found in Reference sohrab\_disloc, . We also note that, for molybdenum, within density functional theory, the hard core configuration was not stable and therefore the stable structure found within MGPT was used as the reference state.

Figure 5 shows the results. Most noticeably, the atomistic landscapes are three times stiffer than the ab initio landscapes. This raises the question whether the approximate factor of three overestimate of theoretical calculations over the extrapolation of the experimental Peierls stresses to zero-temperature is due to defects in the inter-atomic potentials or to failures in the connection between the experiments and the theoretical calculations.

### V.3 Verification of cell size

To compute the Peierls stress, we shall employ a cell of dimensions 23Å$`\times `$ 12Å. To verify that long-range electronic structure effects in metals do not interfere with results in such a cell, we compare the core structure reported previouslysohrab\_disloc  for this cell with a new calculation using a larger cell. Figure 6 shows the result for the core structure in cell of size 41Å$`\times `$ 20Å. The quadrupole array has size $`\stackrel{}{r}_1=9a[1,\overline{1},0],\stackrel{}{r}_2=(5/2)a[1,1,\overline{2}],\text{and}\stackrel{}{r}_3=a[1,1,1]/2`$. The lattice vectors are $`\stackrel{}{a}_1=\stackrel{}{r}_1/2\stackrel{}{r}_2+\stackrel{}{r}_3/2,\stackrel{}{a}_2=\stackrel{}{r}_1/2+\stackrel{}{r}_2+\stackrel{}{r}_3/2\text{ and }\stackrel{}{a}_3=\stackrel{}{r}_3`$. Here we use eight special $`k`$-points kpts2 , which is sufficient as the lattice vectors in the plane of the dislocation have doubled.

We note that the core structure is very similar to previous resultssohrab\_disloc  which used a smaller cell equal in size to that we shall employ for our calculation of the Peierls stress. We therefore do not expect the long-range nature of electronic effects in metallic systems to greatly affect the value which shall extract for the Peierls stress. We also note that the empirical potential results for tantalum (Section IV) had a very large cutoff of 9Å and accurate results were obtained in the smallest cells used in those calculations. These facts lend confidence in the reliability of our density-functional theory predictions below.

### V.4 Results for the Peierls Stress

Figure 7 shows our density functional theory results for energy as a function of strain following the procedure of Section IV.1. We regard each data point as fully relaxed when the magnitude of the residual force on each atom is less than 0.005 eV/Å.

From Figure 7 it is apparent that, at a strain of 0.05 the Peierls stress has been exceeded, whereas at a strain of 0.04 the energy curve lies on the elastic solution. From the curvature of the data in the elastic region ($`C^{}`$ = 35.2 GPa) and these critical strains, we bound the Peierls stress ($`\sigma _{xz}`$) to be between 1.41 GPa and 1.76 GPa. These results are slightly lower than the previous density functional results of $``$ 1.8 GPa obtained using Greens function boundary conditionswoodward2 . However both results are still in reasonable agreement, as we expect our results to be within 20% of the infinite cell limit, while it is unclear how the domain boundary affected the value of the Peierls stress in the Greens function boundary condition calculation.

Figures 8 and 9 show differential displacement mapsvitek  of the dislocation configuration at strains of 0.04 and 0.05, respectively. The solid triangles in the figures represent the location of the center of the dislocation before the application of stain and the upside-down, solid triangle represents the center of the dislocation once it has moved. Figure 8 shows the relaxed dislocation core to remain in the position of its unstressed state, whereas Figure 9 shows that at an applied strain of 0.05, the dislocation moves onto the next triad of columns of atoms. This glide is consistent with such screw dislocations as the initial displacement is along the (01$`\overline{1}`$)-plane. The subsequent motion should be along the (1$`\overline{1}`$0)-plane, and hence overall motion is along a {112}-plane (twinning direction), consistent with the previous density functional theory woodward2  calculations.

Intriguingly, despite the fact that the ab initio energy landscape is less corrugated than that of the inter-atomic potential by a factor of nearly three (Figure 5), the above ab initio result for the Peierls stress is over a factor of two larger than the empirical potential result (Table 1). Moreover, our result is over a factor of five larger than extrapolations of experimental data to zero-temperaturetang , a discrepancy much larger than any effects from the use of either our relatively small quadrupolar cell or the local density approximation to density-functional theory. Certainly, as Section III shows, some of this discrepancy can be due to the extrapolation of experimental data to zero temperature. However, errors in the zero-temperature extrapolation may not likely account for all of the discrepancy as the underestimation in Figure 2 is only $`1020\%`$. One must therefore consider the possibility of other factors such as the effects of kinks, defects, surfaces, strain-rate dependencies and screening from the presence of other dislocations, to accurately predict the experimental results.

## VI Conclusion

This work explores various aspects of the the Peierls stress for the $`111`$ screw dislocation in bcc tantalum. The first non-zero temperature results for the Peierls stress in this system shows both a strong orientation and temperature dependent response, consistent with experimental results. This data also demonstrate that common extrapolations of experimental data tend to underestimate the zero-temperature limit.

We have also presented the first density functional theory calculation for the Peierls stress within perioidic boundary conditions, the approach best suited to metallic systems. The value we find for the Peierls stress is substantially larger than both the experimental extrapolations and current empirical potential results. This difference is much larger than errors which could normally be attributed to the use of a relatively small unit cell or the local density approximation to density-functional theory. The error is also significantly larger than the underestimation we have seen in the the extrapolation of nonzero-temperature data to zero temperature, thus supporting the notion that mechanisms other than the simple Peierls resistance may play an important role in controlling the process of low temperature slip.

## VII Acknowledgments

The authors would like to thank Guofeng Wang for providing us with the new parameters for qEAM potential. This work was supported by an ASCI ASAP Level 2 grant (Contracts No. B338297 and No. B347887). Computational support on ASCI Blue-Pacific provided through the Cal-Tech DOE ASCI center. We thank the members of the H-division at Lawrence Livermore National Laboratories for providing the Ta pseudopotential, the Mo MGPT code, and many useful discussions.


format: text

text:
# Fisher Information With Respect to Cumulants

## I Introduction

Fisher information stats ; fisherinfo  constitutes a central concept in statistical estimation theory which furnishes a variety of useful estimates of deterministic parameters from statistical observations typical of a physical experiment. Its inverse yields a lower bound, called the Cramer-Rao lower bound (CRLB), on the variance of any unbiased estimator of a continuous parameter and thus limits the best precision with which the parameter can be extracted from statistical measurements fisherinfo ; stats .

It is useful to regard a probability density function (pdf) of a random variable as being implicitly parameterized in terms of a translational location parameter, e.g., its mean, median, or mode. The Fisher information relative to such a purely translation parameter is easily seen to be independent of that parameter, and may be defined as the Fisher information of the random variable fisherinfo  itself. The notion of Fisher information of a random variable has been applied to the case in which the random variable is the sample based mean. In the limit of large sample size, asymptotic estimates fishersamplemean  have been obtained in terms of the cumulants of the underlying pdf and their derivatives.

In this correspondence, we present a further generalization of the Fisher information of a continuous random variable. Instead of using only implicit parameterizations, we explicitly parameterize all smooth, well-behaved pdf’s in terms of their cumulants mathworldcumulants ; Kenney1951cumulants . Since the set of all cumulants of a pdf uniquely and completely specifies the pdf, the Fisher information matrix relative to all the cumulants should represent, in effect, the fidelity of estimation of the full pdf from data. The choice of cumulants to parameterize a pdf is a particularly convenient one since, as we shall see, it leads to a simple analytical form for the Fisher information matrix. These concepts can be generalized still further, as we shall argue, to the case of a discrete, integer-valued random variable. We shall present some useful properties of these information measures, discuss their generating functions, and illustrate our considerations with simple examples.

## II Fisher Information of a Continuous Random Variable

Given a continuous random variable $`X`$ distributed according to the pdf $`p(x)>0`$ for all real $`x`$ (where $`x`$ is a statistical realization of $`X`$), we may define a parameterized version of this pdf, $`p(x|\theta )=p(x\theta )`$. Note that $`p(x|\theta )=p(x)`$ if $`\theta =0`$. The Fisher information of $`p(x|\theta )`$ with respect to $`\theta `$ is

$`J_\theta `$ $`=`$ $`\left({\displaystyle \frac{\mathrm{ln}p(x|\theta )}{\theta }}\right)^2`$ (1)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xp(x|\theta )\left({\displaystyle \frac{\mathrm{ln}p(x|\theta )}{\theta }}\right)^2`$
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xp(x\theta )\left({\displaystyle \frac{\mathrm{ln}p(x\theta )}{x}}\right)^2`$
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{p(x)}}\left({\displaystyle \frac{dp}{dx}}\right)^2J(X),`$
where the angled brackets in the first line indicate expectation value with respect to $`p(x|\theta )`$, and $`J(X)`$ is the Fisher information of the random variable $`X`$ fisherinfo . Note that because the pdf is positive for all real values of $`x`$, we can integrate over the infinite interval, and $`J_\theta `$ is manifestly independent of $`\theta `$. Thus, any location parameter may be used for $`\theta `$, and $`J(X)`$ is therefore a functional only of the shape of the distribution, independent of its absolute location on any axis.

The Fisher information of a random variable $`J(X)`$ has two well-known interpretations stats ; fisherinfo . First, $`J(X)`$ quantifies the statistical precision with which the location parameter $`\theta `$ of a pdf on which the pdf depends translationally can be estimated from data drawn according to the pdf $`p(x\theta )`$. On the other hand, because $`J(X)`$ measures the mean squared slope of the log-likelihood function $`\mathrm{ln}p(x)`$, it typically correlates with the narrowness of the pdf or, equivalently, with the degree of statistical reproducibility of the values assumed by the variable $`X`$. In the Bayesian context, this narrowness is related to the extent of prior knowledge about $`X`$. These two interpretations are related in that a narrower pdf will provide higher “resolution” when used as a measurement tool for determining the location parameter $`\theta `$.

A question of central interest in this correspondence is the following: What information measures characterize the fidelity of a statistical determination of the full pdf, not just its location parameter? We shall see presently that a particularly simple answer to this question can be obtained in terms of the Fisher information matrix elements relative to the cumulants of the pdf of the random variable.

## III Cumulants of a PDF and Associated Fisher Information

Every pdf $`p(x)`$ has a unique characteristic function associated with it, given by its Fourier transform,

$`M(\nu )e^{i\nu x}={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xp(x)e^{i\nu x},`$ (2)
which, in most cases, may be expressed in a power series in $`\nu `$ in terms of the moments of the pdf:

$`M(\nu )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(i\nu )^n}{n!}}\mu _n^{},`$ (3)
where $`\mu _n^{}x^n`$ is the $`n`$th moment of the pdf (about 0). Writing the logarithm of the characteristic function in a similar series form defines the cumulants $`\kappa _n`$ of the pdf as coefficients in the series expansion mathworldcumulants ; Kenney1951cumulants

$`L(\nu )\mathrm{ln}M(\nu ){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(i\nu )^n}{n!}}\kappa _n.`$ (4)
Note that the function $`L(\nu )`$ may be regarded as the generating function for the cumulants, since

$`\kappa _n={\displaystyle \frac{1}{i^n}}{\displaystyle \frac{d^nL}{d\nu ^n}}|_{\nu =0}.`$ (5)
The cumulants $`\kappa _n`$ are related to the mean and central moments of the pdf, the first few of them taking the form $`\kappa _1=\mu _1^{},\kappa _2=\mu _2,\kappa _3=\mu _3,\kappa _4=\mu _42\mu _2^2,\mathrm{}`$, where $`\mu _1^{}`$ and $`\mu _2`$ are the mean and variance of the pdf, while its third and fourth central moments, $`\mu _3,\mu _4`$, are related to its skewness and kurtosis. Since a pdf is related to its characteristic function by a Fourier transform, we may thus parameterize $`p(x)`$ in terms of the entire collection of its cumulants (indicated by $`𝜿`$) as the following inverse Fourier transform:

$`p(x)p(x|𝜿){\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\nu e^{i\nu x}\mathrm{exp}\left\{{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(i\nu )^n}{n!}}\kappa _n\right\}.`$ (6)
The Fisher information with respect to a set of real estimation parameters $`𝜽\{\theta _1,\theta _2,\mathrm{}\}`$ is defined as a positive semi-definite matrix with elements stats ; fisherinfo

$`J_{mn}^{(𝜽)}`$ $``$ $`{\displaystyle \frac{\mathrm{ln}p(x|𝜽)}{\theta _m}}{\displaystyle \frac{\mathrm{ln}p(x|𝜽)}{\theta _n}},`$ (7)
where the expectation value is taken over the pdf, $`p(x|𝜽)`$, of the statistical data from which the parameters are estimated. The superscript $`(𝜽)`$ indicates the parameter set being used, and both $`n`$ and $`m`$ range over the indices of the parameters in $`𝜽`$. For a continuous pdf $`p(x)`$ that does not vanish at any finite value of $`x`$, this becomes

$`J_{mn}^{(𝜽)}={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{p(x|𝜽)}}{\displaystyle \frac{p(x|𝜽)}{\theta _m}}{\displaystyle \frac{p(x|𝜽)}{\theta _n}}.`$ (8)
Noting from Eq. (6) that

$`{\displaystyle \frac{p(x|𝜿)}{\kappa _n}}={\displaystyle \frac{(1)^n}{n!}}{\displaystyle \frac{^np(x|𝜿)}{x^n}}`$ (9)
and using Eq. (8) with $`𝜽`$ chosen as the cumulant vector $`𝜿`$, we obtain the Fisher information matrix relative to the cumulants,

$`J_{mn}^{(𝜿)}={\displaystyle \frac{(1)^{m+n}}{m!n!}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{p(x|𝜿)}}{\displaystyle \frac{^mp(x|𝜿)}{x^m}}{\displaystyle \frac{^np(x|𝜿)}{x^n}},m,n=0,1,2,\mathrm{}.`$ (10)
While $`\kappa _0`$ is not a true cumulant, including the possibility that $`n`$ and/or $`m=0`$ in Eq. (10) leads formally to a particularly convenient generating function for these elements, as we shall see in Sec. IV.

Equation (10) is the most important result of this correspondence. It defines a complete hierarchy matrix of information measures, which we may call the cumulant information matrix (CIM). The diagonal elements of the inverse of the CIM yield the full hierarchy of CRLB’s on how precisely the various cumulants of $`p(x|𝜿)`$ may be estimated from the statistical realizations of $`X`$ stats  and thus, in a sense, how well the entirety of the pdf may be estimated.

For $`m=n=1`$, the CIM element $`J_{11}^{(𝜿)}`$ reduces to $`J(X)`$ defined in Eq. (1) as the Fisher information of the random variable. This result confirms our earlier interpretation of $`J(X)`$ as the precision with which a fiducial location parameter of the pdf may be estimated. Without knowledge of any higher order cumulants, it is the first cumulant (the mean) that furnishes the most useful location parameter of a pdf.

## IV A Generating Function for the Cumulant Information Matrix

A useful technique for evaluating the CIM elements is the method of generating functions. Upon multiplying both sides of Eq. (10) by $`\lambda ^m\mu ^n`$, then summing over all non-negative integral values of $`m`$ and $`n`$, and finally interchanging the order of integration and summations, a procedure justified by the uniform convergence of the involved Taylor expansions, we obtain the following result:

$`J^{(𝜿)}(\lambda ,\mu )`$ $``$ $`{\displaystyle \underset{m,n=0}{\overset{\mathrm{}}{}}}J_{mn}^{(𝜿)}\lambda ^m\mu ^n`$ (11)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{p(x|𝜿)}}\left[{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\lambda )^m}{m!}}{\displaystyle \frac{^mp(x|𝜿)}{x^m}}\right]\left[{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(\mu )^n}{n!}}{\displaystyle \frac{^np(x|𝜿)}{x^n}}\right].`$ (12)
The two pairs of square brackets in Eq. (12) enclose the Taylor expansions of $`p(x\lambda |𝜿)`$ and $`p(x\mu |𝜿)`$, respectively. We thus arrive at a rather simple form of the function $`J^{(𝜿)}(\lambda ,\mu )`$:

$`J^{(𝜿)}(\lambda ,\mu )={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle \frac{p(x\lambda |𝜿)p(x\mu |𝜿)}{p(x|𝜿)}}.`$ (13)
The function $`J^{(𝜿)}(\lambda ,\mu )`$ given by Eq. (13) is a generating function for the CIM elements, since from Eq. (11) an arbitrary element $`J_{mn}^{(𝜿)}`$ may be expressed as its partial derivative

$`J_{mn}^{(𝜿)}={\displaystyle \frac{1}{m!n!}}{\displaystyle \frac{^{m+n}}{\lambda ^m\mu ^n}}J^{(𝜿)}(\lambda ,\mu )|_{\lambda =\mu =0}.`$ (14)
The elements $`J_{0n}^{(𝜿)}=J_{n0}^{(𝜿)}`$, as mentioned earlier, do not correspond to any information relative to any cumulants, and can, in fact, be easily shown to vanish for $`n1`$, while $`J_{00}^{(𝜿)}=1`$, the pdf normalization.

## V Example: Gaussian pdf

The Gaussian pdf provides an analytically tractable illustration of the results of this correspondence. For the Gaussian pdf,

$`p(x)={\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}e^{(xx_0)^2/2\sigma ^2},`$ (15)
its characteristic function is also Gaussian, and its cumulant generating function is thus quadratic,

$`L(\nu )=\mathrm{ln}M(\nu )=i\nu x_0\nu ^2\sigma ^2/2.`$ (16)
The first two cumulants are thus its mean $`x_0`$ and variance $`\sigma ^2`$, while all higher order cumulants vanish identically. Parameterizing the Gaussian in terms of its cumulants $`𝜿`$, we obtain, very simply,

$`p(x|𝜿)={\displaystyle \frac{1}{\sqrt{2\pi \kappa _2}}}e^{(x\kappa _1)^2/2\kappa _2}.`$ (17)
The CIM generating function can be easily evaluated for the Gaussian pdf (17), for which the integrand in Eq. (13) is also Gaussian and easily integrated, with the result

$`J^{(𝜿)}(\lambda ,\mu )=e^{\lambda \mu /\kappa _2}=e^{\lambda \mu /\sigma ^2}.`$ (18)
The individual CIM elements then follow from a use of Eq. (14):

$`J_{mn}^{(𝜿)}={\displaystyle \frac{\delta _{mn}}{n!\sigma ^{2n}}}.`$ (19)
The diagonal nature of the CIM implies simple CRLB’s on the variance of any set of unbiased estimators, $`\{\widehat{\kappa }_n\}`$, of the cumulants of the pdf:

$`\text{var}(\widehat{\kappa }_n)(\widehat{\kappa }_n\kappa _n)^2n!\sigma ^{2n}.`$ (20)
For general, biased estimators, the right-hand side of Eq. (20) must be replaced by the $`nn`$diagonal element of the matrix $`B^TIB`$, where the matrix $`I`$ is the inverse of the CIM, which for the Gaussian case is diagonal with the $`nn`$element equal to $`n!\sigma ^{2n}`$, $`B`$ is the bias matrix, with elements $`B_{mn}=\widehat{\kappa }_m/\kappa _n`$, and $`B^T`$ is its transpose stats .

The $`k`$-statistics furnish useful unbiased estimators of the cumulants of a distribution based on a finite sample drawn from that distribution mathworldk-stats ; Kenney1962 ; Halmos1946 . It is well known that no cumulant estimator exists with a smaller variance than that of the corresponding $`k`$-statistic. What our results show is that the $`k`$-statistics for a Gaussian pdf are also asymptotically efficient estimators in the sense that the CRLB (20) is achieved in the limit that the sample size $`N\mathrm{}`$. This asymptotic efficiency has not been previously derived for any of the $`k`$-statistics of order higher than 2. The variances of the first few $`k`$-statistics for the Gaussian pdf are given below, along with the associated CRLB’s from the CIM analysis, the latter given as the rightmost terms in the following expressions:

$`\text{var}(k_1)`$ $`=`$ $`{\displaystyle \frac{\kappa _2}{N}}={\displaystyle \frac{\sigma ^2}{N}}`$ (21a)
$`\text{var}(k_2)`$ $`=`$ $`{\displaystyle \frac{2\kappa _2^2}{N1}}{\displaystyle \frac{2\sigma ^4}{N}}`$ (21b)
$`\text{var}(k_3)`$ $`=`$ $`{\displaystyle \frac{6N\kappa _2^3}{(N1)(N2)}}{\displaystyle \frac{6\sigma ^6}{N}}`$ (21c)
$`\text{var}(k_4)`$ $`=`$ $`{\displaystyle \frac{24(N+1)N\kappa _2^4}{(N1)(N2)(N3)}}{\displaystyle \frac{24\sigma ^8}{N}}`$
$`\mathrm{}`$
(The CRLB for $`N`$ trials is $`1/N`$ times the CRLB for one trial.) Notice that for the Gaussian, $`k_1`$ is an efficient estimator for any $`N`$, and the others are asymptotically efficient (indicated by $``$) in the limit $`N\mathrm{}`$. This is true of all higher-order $`k`$-statistics for the Gaussian, as well.

Two remarks are in order here. First, although any unbiased estimator of a third or higher order cumulant is, on average, zero for the Gaussian pdf, there is a finite statistical scatter in the data from which the cumulants, regardless of their order, are estimated. Second, the sharp increase of the CRLB’s (20) with increasing $`n`$ is a reflection of the sharply decreased probability of occurrence of values of a Gaussian variate in the wings of the pdf to which the higher order cumulants are increasingly sensitive as a function of their order. This is in fact a general result for any localized distribution; the CIM elements $`0`$ as $`n,m\mathrm{}`$, because of the factor of $`n!m!`$ in the denominator of Eq. (10), resulting in CRLB’s that increase without bound as the order of the cumulant being estimated increases.

## VI An Inverse CIM Generating Function

The CIM elements are of value only insofar as they give a general sense of how much information the data contain about the cumulants of the pdf from which the data are drawn. It is the inverse of the CIM that is needed to establish the CRLB, and such an inverse is often difficult to calculate. If a less stringent bound is acceptable, one can use the reciprocal of the $`nn`$-diagonal element of the CIM diagonalbound  to bound the minimum variance from below. Such a bound, though easy to write down, is not optimal, however, since the CRLB represents in general a greater lower bound. The desire for a simple way to calculate the inverse CIM motivates the following attempt to define an appropriate generating function for the inverse CIM and relate it to the CIM generating function defined in Eq. (13).

Let us define the generating function $`I^{(𝜿)}(\nu ,\eta )`$ for the matrix elements $`I_{mn}^{(𝜿)}`$ of the inverse of the CIM as

$`I^{(𝜿)}(\nu ,\eta ){\displaystyle \underset{m,n=0}{\overset{\mathrm{}}{}}}I_{mn}^{(𝜿)}{\displaystyle \frac{\nu ^m\eta ^n}{m!n!}},`$ (22)
from which the inverse CIM elements may be obtained as follows:

$`I_{mn}^{(𝜿)}={\displaystyle \frac{^{m+n}}{\nu ^m\eta ^n}}I^{(𝜿)}(\nu ,\eta )|_{\nu =\eta =0}.`$ (23)
The factorials, $`m!`$ and $`n!`$, are included in the definition (22) to allow appropriate convergence of the generating function (Compare Eq. (11)). We may also define two “marginal” generating functions,

$`J_n^{(𝜿)}(\lambda ){\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}J_{mn}^{(𝜿)}\lambda ^m\text{and}I_n^{(𝜿)}(\nu ){\displaystyle \underset{\mathrm{}=0}{\overset{\mathrm{}}{}}}I_\mathrm{}n^{(𝜿)}{\displaystyle \frac{\nu ^{\mathrm{}}}{\mathrm{}!}},`$ (24)
that are related to the full generating functions (13) and (22) by

$`J_n^{(𝜿)}(\lambda )={\displaystyle \frac{1}{n!}}{\displaystyle \frac{^n}{\mu ^n}}J^{(𝜿)}(\lambda ,\mu )|_{\mu =0}\text{and}I_n^{(𝜿)}(\nu )={\displaystyle \frac{^n}{\eta ^n}}I^{(𝜿)}(\nu ,\eta )|_{\eta =0}.`$ (25)
By multiplying the two marginals (24), summing over the index $`n`$ from 0 to $`\mathrm{}`$, and noting that the symmetric matrices $`I_{mn}^{(𝜿)}`$ and $`J_{mn}^{(𝜿)}`$ are mutual inverses, we obtain the following relation:

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}J_n^{(𝜿)}(\lambda )I_n^{(𝜿)}(\nu )`$ $`=`$ $`{\displaystyle \underset{k,m=0}{\overset{\mathrm{}}{}}}\lambda ^k{\displaystyle \frac{\nu ^m}{m!}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}J_{kn}^{(𝜿)}I_{mn}^{(𝜿)}`$ (26)
$`=`$ $`{\displaystyle \underset{k,m=0}{\overset{\mathrm{}}{}}}\lambda ^k{\displaystyle \frac{\nu ^m}{m!}}\delta _{km}`$
$`=`$ $`e^{\lambda \nu }.`$
A Fourier analysis gives a second, integral relation,

$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\mu {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\eta e^{i\mu \eta }J^{(𝜿)}(\lambda ,i\mu )I^{(𝜿)}(\nu ,\eta )=e^{\lambda \nu },`$ (27)
as shown in the Appendix. These relations are valid for all values of $`\lambda `$ and $`\nu `$ in the complex plane. In general, relations (26) and (27) cannot be solved analytically. The case of the Gaussian pdf is an exception, for which the preceding relations can be solved and the generating function $`I^{(𝜿)}(\nu ,\eta )`$ derived exactly,

$`I^{(𝜿)}(\nu ,\eta )=e^{\nu \eta \sigma ^2}.`$ (28)
Applying Eq. (23) to this result yields the same CRLB’s (20) as obtained in Sec. V. While the CIM for the Gaussian is admittedly trivial to invert, and thus it is not necessary to use the generating function method in this case, we anticipate that in some cases it will prove easier to approximate the solution of either relation (26) or (27) than it would be to invert the full CIM.

## VII Generalization to a Discrete Random Variable

While the Fisher information of a continuously random variable over the infinite interval may be defined through the use of an arbitrary location parameterization, discrete distributions do not readily permit an analogous definition. Although one cannot treat a discrete random variable as a continuously-valued parameter, the class of discrete distributions defined positive over the entirety of the integers can always be parameterized in terms of their cumulants and a corresponding CIM matrix defined.

For a discrete random variable $`X`$ distributed according to $`P(x)`$ with $`x`$ ranging over all integers, we may parameterize this distribution as

$`P(x)P(x|𝜿){\displaystyle \frac{1}{2\pi }}{\displaystyle _\pi ^\pi }𝑑\nu e^{i\nu x}\mathrm{exp}\left\{{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(i\nu )^n}{n!}}\kappa _n\right\},`$ (29)
where

$`\kappa _n{\displaystyle \frac{1}{i^n}}{\displaystyle \frac{d^n\mathrm{ln}e^{i\nu x}}{d\nu ^n}}|_{\nu =0}`$ (30)
are the cumulants in exact analogy to the continuous case. The only formal difference is in terms of the limits of integration in (29), reflecting the fact that $`e^{i\nu x}`$ is now a Fourier series.<sup>1</sup><sup>1</sup>1 For a random variable that has a one-sided discrete realization over the set of all non-negative integers only, the characteristic function $`\mathrm{exp}(i\nu x)`$ is the $`z`$-transform Bracewell99  of the distribution $`P(x)`$, rather than its Fourier series, with $`z=\mathrm{exp}(i\nu )`$. Notwithstanding this difference, the parameterization of the distribution in terms of its cumulants is formally identical to that given in Eq. (29), and the considerations of this section apply essentially unchanged.

The parameterized distribution $`P(x|𝜿)`$ coincides with $`P(x)`$ for all integral values of $`x`$, but it is also $`C^{\mathrm{}}`$ and therefore defined for intermediate values of $`x`$ as well. These intermediate values are not probabilities of anything, and they can range outside the interval $`[0,1]`$. Still, because of its coincidence with $`P(x)`$ at integral $`x`$, $`P(x|𝜿)`$ may be substituted for $`P(x)`$ when calculating an expectation value of any function of $`X`$. Most importantly, though, we may now also take derivatives of this parameterized distribution with respect to $`x`$, which are related to partial derivatives with respect to the distribution’s cumulants exactly as in the continuous case, i.e.

$`{\displaystyle \frac{P(x|𝜿)}{\kappa _n}}={\displaystyle \frac{(1)^n}{n!}}{\displaystyle \frac{^nP(x|𝜿)}{x^n}}.`$ (31)
This allows us to define the discrete CIM elements exactly as in the continuous case, with the infinite integral replaced by an infinite sum:

$`J_{mn}^{(𝜿)}={\displaystyle \frac{(1)^{m+n}}{m!n!}}{\displaystyle \underset{x=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{P(x|𝜿)}}{\displaystyle \frac{^mP(x|𝜿)}{x^m}}{\displaystyle \frac{^nP(x|𝜿)}{x^n}},m,n=0,1,2,\mathrm{}.`$ (32)
The possibility that $`n`$ or $`m=0`$ is also included in this definition for assistance in defining a generating function for these matrix elements, namely

$`J^{(𝜿)}(\lambda ,\mu )={\displaystyle \underset{x=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \frac{P(x\lambda |𝜿)P(x\mu |𝜿)}{P(x|𝜿)}},`$ (33)
in exact analogy to the continuous case. Without modification, Eq. (14) may be used to generate the individual CIM elements, and the generating function relations given in Eqs. (26) and (27) also hold.

This is a powerful result. While it is not possible to directly define the Fisher information of a discrete random variable as in the continuous case, we can define an analogous quantity: the Fisher information with respect to the mean of the cumulant-parameterized distribution. This is $`J_{11}^{(𝜿)}`$ of the discrete CIM matrix. Since in the continuous case $`J_{11}^{(𝜿)}=J(X)`$, we can define $`J(X)J_{11}^{(𝜿)}`$ for a discrete distribution. The CIM matrix provides a complete hierarchy of information measures with respect to the cumulants in the discrete case, as well, and the inverse of the CIM again gives the CRLB on any unbiased estimators of the cumulants.

## VIII Concluding Remarks

The Fisher information of a continuous random variable can be interpreted as the fidelity with which a fiducial location parameter of a pdf (such as the mean) may be estimated from statistical data drawn according to that distribution. In this correspondence, we have introduced a more robust information measure—the cumulant information matrix (CIM)—whose inverse bounds the variance of any estimates of the cumulants of a pdf and, consequently, the fidelity with which the entire pdf may be estimated. The Fisher information of the random variable is included in this measure. We have also extended the CIM concept to discrete random variables defined over integers, for which the notion of Fisher information of the random variable is ill defined. We have also derived a generating function for the CIM and given two relations between this generating function and a generating function for the inverse CIM, which we hope will prove useful in calculating CRLB’s. Further generalizations of this work could include defining a cumulant information matrix or an analogous quantity for multivariate distributions. (See McCullagh McCullagh  for information on multivariate cumulants.)

We noted in Sections III and VII that the probability distribution to be parameterized must be strictly positive. This was done to avoid singularities in the Fisher information calculations. As a general rule, the CIM for an arbitrary random variable must be defined by restricting the interval of integration or summation to the actual sample space of that variable. This is particularly important when a priori constraints like finiteness of support place restrictions on the set of possible values of the allowed cumulant vectors. One constraint that is implicitly active in all of our calculations is that the probability distribution, whether continuous or discrete, is restricted to nonnegative values. While this is obvious from the standpoint of probability theory, it must be explicitly imposed on the space of all estimable probability distributions. One only need consider the parameterization (6) to see why the cumulants $`𝜿`$ cannot take arbitrary values if the left hand side of that equation is to remain nonnegative. CRLB’s calculated in the present paper apply only to estimators that do not explore the boundaries of the parameter space beyond which negative distributions are encountered. Incorporating edge effects will, in general, reduce the minimum variance of a more general estimator, but such bounds are not calculable from the Fisher information matrix since inequality constraints do not affect its form. However, certain other constraints such as support constraints can alter the Fisher information matrix and change the method that one must use to calculate CRLB’s. The effects of constraints on CRLB’s are explored in detail by Gorman and Hero Gorman1990 .

The considerations of this paper may be relevant to the general area of inverse problems such as those concerning image restoration from noisy image data. A particularly useful viewpoint to adopt in discussions of image processing is to treat the spatial distribution of intensity in an image, when properly normalized, as representing the probability distribution of the emission or detection of a photon over the image. From this perspective, image restoration is equivalent to the problem of estimating a probability distribution, the very problem we have discussed here. The presence of noise in the actual image data greatly compounds this estimation problem, a subject that requires further study. A noteworthy algorithm which makes essential use of this statistical viewpoint is the maximum entropy method Narayan86 . Blind deconvolution methods Schulz93  provide another example. They rely on the existence of constraints to recover the point-spread function as well as the source intensity distribution, both of which may be regarded as appropriate pdf’s.

## Acknowledgments

The work reported here was supported in part by the US Air Force Office of Scientific Research under grants F49620-00-1-0155 and F49620-01-1-0321. We are also pleased to receive from Dr. Chris Lloyd a reprint of his paper cited here.

## Appendix A Derivation of A Generating Function Relation

To prove Eq. (27), we first note that

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ye^{ixy}x^ny^m`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xx^n{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑yy^me^{ixy}`$ (34)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xx^ni^m{\displaystyle \frac{^m}{x^m}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ye^{ixy}`$
$`=`$ $`2\pi i^m{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑xx^n\delta ^{(m)}(x),`$
where $`\delta ^{(m)}(x)`$ is the $`m`$th derivative of the Dirac $`\delta `$-function. Integrating by parts $`m`$ times in the right-hand side of Eq. (34) reduces that integral to $`(1)^m`$ times the $`m`$th derivative of $`x^n`$ evaluated at $`x=0`$, namely to $`(1)^mm!\delta _{mn}`$, and the following identity results:

$`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑ye^{ixy}{\displaystyle \frac{(ix)^n}{n!}}y^m`$ $`=`$ $`\delta _{mn}.`$ (35)
Now, from the generating function definitions (11) and (22) and from the marginal definitions (24), it is clear that

$`J^{(𝜿)}(\lambda ,\mu )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}J_n^{(𝜿)}(\lambda )\mu ^n\text{and}I^{(𝜿)}(\nu ,\eta )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}I_n^{(𝜿)}(\nu ){\displaystyle \frac{\eta ^n}{n!}}.`$ (36)
Plugging these results into the left-hand side (LHS) of Eq. (27), we obtain

LHS $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\mu {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\eta e^{i\mu \eta }\left[{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}J_m^{(𝜿)}(\lambda )(i\mu )^m\right]\left[{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}I_n^{(𝜿)}(\nu ){\displaystyle \frac{\eta ^n}{n!}}\right]`$ (37)
$`=`$ $`{\displaystyle \underset{m,n=0}{\overset{\mathrm{}}{}}}J_m^{(𝜿)}(\lambda )I_n^{(𝜿)}(\nu )\left[{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\mu {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\eta e^{i\mu \eta }(i\mu )^m{\displaystyle \frac{\eta ^n}{n!}}\right]`$
$`\stackrel{(a)}{=}`$ $`{\displaystyle \underset{m,n=0}{\overset{\mathrm{}}{}}}J_m^{(𝜿)}(\lambda )I_n^{(𝜿)}(\nu )\delta _{mn}`$
$`\stackrel{(b)}{=}`$ $`e^{\lambda \nu }`$
where Eq. (35) was used in step (a), and Eq. (26) allowed step (b). Thus, relation (27) is established.