text:
# Fluctuation-dissipation theorem and flux noise in overdamped Josephson junction arrays

## I Introduction

The two basic experimental methods used for contactless investigation of finite frequency properties of two-dimensional superconducting systems (such as thin films , Josephson junction arrays  and wire networks ) are two-coil mutual-inductance technique  and flux noise power spectrum analysis . The first of them is based on measurement of the voltage induced in the detection coil by the currents flowing in the sample under the action of ac electric field produced by the current in the other (driving) coil. For the given geometry of the coils the measured signal can be used to extract  the complex frequency dependent sheet impedance $`Z_{\mathrm{}}(\omega )`$ of the sample on the assumption that for wave-lengths larger than the characteristic dimensions of the detection coil $`Z_{\mathrm{}}(\omega )`$ is not wave-length dependent.

In the case of the flux noise spectrum analysis the approaches to interpretation of the experimental data are much more varied. The theoretical predictions of the flux noise spectrum used for comparison with experimental data are found by relating it with  or (no less often) by replacing it by  a correlation function describing the vortex distribution and, naturally, turn out to be dependent on the particular choice of assumptions concerning the form of this distribution. Numerical simulations also demonstrate a clear tendency towards studying the vortex number noise  rather then the flux noise. The only attempt to achieve a description of the flux noise power spectrum in terms of the sheet impedance of the sample taking into account the actual geometry of the detection coil has been undertaken by Kim and Minnhagen . However, this calculation is also based on expressing all quantities in terms of the vortex gas correlation functions and, therefore, can be trusted only in a limited range of parameters.

In the present work we argue that in the case of a resistively shunted Josephson junction array the general expression for the flux noise spectrum can be found without artificial decomposition of all fluctuations into the vortex part (which is usually assumed to be responsible for the flux noise) and the remaining so-called ”spin-wave” part (which is traditionally neglected). Although in semi-phenomenological treatment  of two-dimensional superconductors such decomposition seems to be inevitable, the case of an overdamped Josephson junction array allows for application of a more universal approach to calculation of the flux noise power spectrum. It is based on the direct relation (discussed in Sec. 2) of the flux noise with current correlations, which, on the other hand, can be expressed in terms of the complex frequency dependent sheet impedance with the help of the fluctuation-dissipation theorem.

The additional advantage of such approach is that it allows to include into consideration in systematic way the mutual influence between magnetic field fluctuations and current fluctuations (the screening effects), which insofar has been neglected in theoretical works  devoted to flux noise spectrum analysis. The form of the Hamiltonian, which should be used for the description of a resistively shunted array in presence of self-induced magnetic fields, is discussed in Sec. 3, and the corresponding dynamic equations in Sec. 4.

The explicit form of the fluctuation-dissipation theorem for resistively shunted Josephson junction array is derived in Sec. 5. It shows that the current correlations in the array are determined by the response of the current to the external electric field and not directly by the sheet impedance of the array (which is defined as a response to the total electric field) The nature of the expression for currents correlation function, related with the peculiarities of the two-dimensional geometry, allows to expect the same expression to be applicable for arbitrary two-dimensional systems in which capacitive effects can be neglected.

Our main result, the relation between the flux noise power spectrum and the frequency dependent sheet impedance of a two-dimensional superconductor, is presented in Sec. 6, which includes also the discusssion of the noise spectrum dependence on the parameters of the detection coil and comparison of the results with those of other authors.

## II Flux noise and current correlations

In a flux noise experiment one measures and analyses the time dependence of a voltage created in a detection coil by fluctuations of currents in some conducting (or superconducting) object. This voltage is determined by the time derivative of the magnetic flux penetrating the coil, and the value of the flux can be expressed in terms of the current density distribution $`𝐣(𝐫)`$ \[$`𝐫=(x_1,x_2,x_3)`$\] inside the object with the help of the Biot-Savart’s law, which in the Coulomb gauge ($`\text{div}𝐀=0`$) can be written as

$$\mathrm{\Delta }𝐀(𝐫)=\mu _0𝐣^t(𝐫),$$
(1)
where $`𝐀(𝐫)`$ is the vector potential defining the distribution of the magnetic field (magnetic induction) $`𝐁(𝐫)=\text{rot}𝐀(𝐫)`$ created by $`𝐣(𝐫)`$,

$$\mathrm{\Delta }\frac{^2}{x_1^2}+\frac{^2}{x_2^2}+\frac{^2}{x_3^2}$$
(2)
is the three-dimensional Laplacian and

$$𝐣^t(𝐫)𝐣(𝐫)\mathrm{\Delta }^1\text{grad}\text{div}𝐣(𝐫)$$
(3)
is the transverse part of $`𝐣(𝐫)`$. Magnetic fields produced by the longitudinal part of $`𝐣(𝐫)`$ cancel each other.

In the case of a system which can be considered as effectively two-dimensional and situated (for simplicity) in the plane $`x_3=0`$ the three-dimensional current density $`𝐣(𝐫)`$ is reduced to

$$𝐣(𝐫)=𝐢(𝐱)\delta (x_3),$$
(4)
where $`𝐱x_\alpha `$ ($`\alpha =1,2`$) is the two-dimensional vector defining the position of a point in the plane $`x_3=0`$ and $`𝐢i_\alpha `$ is the two-dimensional vector describing the two-dimensional current density.

Substitution of Eq. (4) into Eq. (1) allows then to find that

$`𝐀(𝐫)`$ $`=`$ $`\mu _0{\displaystyle \frac{d^2𝐪}{(2\pi )^2}\frac{dq_3}{2\pi }\frac{\mathrm{exp}i(\mathrm{𝐪𝐱}+q_3x_3)}{q^2+q_3^2}𝐢^t(𝐪)}`$ (5)
$`=`$ $`{\displaystyle \frac{\mu _0}{2}}{\displaystyle \frac{d^2𝐪}{(2\pi )^2}\frac{\mathrm{exp}(i\mathrm{𝐪𝐱}q|x_3|)}{q}𝐢^t(𝐪)},`$ (6)
where

$$𝐢(𝐪)=d^2𝐱\mathrm{exp}(i\mathrm{𝐪𝐱})𝐢(𝐱)$$
(7)
is the (two-dimensional) Fourier transform of $`𝐢(𝐱)`$, $`q=|𝐪|`$ and

$$𝐢^t(𝐪)𝐢(𝐪)\widehat{𝐪}(\widehat{𝐪}𝐢(𝐪))$$
(8)
is the transverse part of $`𝐢(𝐪)`$, $`\widehat{𝐪}𝐪/q`$ being the unit vector parallel to $`𝐪`$.

In the simplest case a coil can be approximated by a closed circular ring. Integration of $`𝐀(𝐫)`$ over the perimeter of the ring $`x_1^2+x_2^2=r^2`$ situated at the distance $`h`$ from the plane $`x_3=0`$ gives

$$\mathrm{\Phi }=𝑑𝐫𝐀(𝐫)=\frac{\mu _0}{2}\frac{d^2𝐪}{(2\pi )^2}F(𝐪)i^t(𝐪),$$
(9)
where

$$i^t(𝐪)=\underset{\alpha ,\beta }{}ϵ_{\alpha \beta }\widehat{q}_\alpha i_\beta (𝐪)$$
(10)
is the amplitude of $`𝐢^t(𝐪)`$, $`ϵ_{\alpha \beta }`$ is the unit antisymmetric tensor,

$$F(𝐪)=\frac{2\pi rJ_1(qr)}{q}\mathrm{exp}(qh)$$
(11)
is the geometrical factor depending on the parameters of the coil and $`J_1(z)`$ is the first order Bessel function. In the case when the coil can be considered as consisting of $`N`$ turns separated by the distance $`b`$ from each other, the expression for $`F(𝐪)`$ should also include an additional factor obtained by summation of contributions from different turns :

$$F(𝐪)=\frac{2\pi rJ_1(qr)}{q}\mathrm{exp}(qh)\frac{1\mathrm{exp}(Nqb)}{1\mathrm{exp}(qb)}.$$
(12)
The power spectrum of the flux noise is given by the flux-flux correlation function

$$S(\omega )=𝑑t\mathrm{\Phi }(t_0+t)\mathrm{\Phi }(t_0)\mathrm{exp}(i\omega t)$$
(13)
and with the help of Eq. (9) can be expressed in terms of the current density correlation function

$`𝐢^t𝐢^t_{𝐪\omega }{\displaystyle d^2𝐱}`$ $`dt`$ $`[\mathrm{exp}(i\mathrm{𝐪𝐱}+i\omega t)`$ (14)
$`\times `$ $`𝐢^t(𝐱_\mathrm{𝟎}+𝐱,t_0+t)𝐢^t(𝐱_\mathrm{𝟎},t_0)]`$ (15)
as

$$S(\omega )=\frac{\mu _0^2}{4}\frac{d^2𝐪}{(2\pi )^2}F^2(𝐪)𝐢^t𝐢^t_{𝐪\omega }.$$
(16)
## III The Hamiltonian of a Josephson junction array

When self-induced magnetic field is taken into account a square Josephson junction array can be described by the Hamiltonian

$`H`$ $`=`$ $`J{\displaystyle \underset{𝐧,\alpha }{}}\mathrm{cos}\left(_\alpha \phi _𝐧A_{𝐧\alpha }\right)`$ (18)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐧,𝐤}{}}(\times A)_𝐧M_{\mathrm{𝐧𝐤}}^1(\times A)_𝐤,`$
where $`\phi _𝐧`$ is the phase of the order parameter on the $`𝐧`$-th superconducting island, the variables $`A_{𝐧\alpha }`$ (defined on the bonds of the lattice) are determined by the integral of the vector potential $`𝐀(𝐫)`$ over the line connecting the geometrical centers of the neighboring superconducting islands:

$$A_{𝐧\alpha }=\frac{2e}{\mathrm{}}_{a𝐧}^{a(𝐧+𝐞_\alpha )}𝑑𝐫𝐀(𝐫)$$
(19)
and $`a`$ is period of the lattice.

The first term in Eq. (18) describes the Josephson energy of the junctions in the array. The coupling constant $`J`$ entering this term is determined by the critical current $`I_c`$ of a single junction:

$$J=\frac{\mathrm{}}{2e}I_c,$$
(20)
which is assumed to be the same for all junctions, whereas $`_\alpha \phi _𝐧`$ denotes the difference of $`\phi _𝐧`$ between the neighboring sites of the lattice:

$$_\alpha \phi _𝐧\phi _{𝐧+𝐞_\alpha }\phi _𝐧.$$
(21)
Here $`𝐞_1=(1,0)`$ and $`𝐞_2=(0,1)`$. Notice that the combination

$$\theta _{𝐧\alpha }_\alpha \phi _𝐧A_{𝐧\alpha }$$
(22)
which enters as the argument of the Josephson energy $`E_J(\theta )=J\mathrm{cos}\theta `$ is a gauge-invariant quantity.

The second term in Eq. (18) is the energy of the magnetic field

$$E_{\mathrm{mf}}=\frac{1}{2\mu _0}d^3𝐫𝐁^2(𝐫)$$
(23)
expressed in terms of the variables $`A_{𝐧\alpha }`$. The matrix $`M_{\mathrm{𝐧𝐤}}M(𝐧𝐤)`$ is usually called the mutual inductance matrix  and

$$(\times A)_𝐧\underset{\alpha ,\beta }{}ϵ_{\beta \alpha }_\beta A_{𝐧\alpha }$$
(24)
is the directed sum of the variables $`A_{𝐣\alpha }`$ along the perimeter of a lattice plaquette (the lattice equivalent of $`\text{rot}𝐀`$) and is proportional to the magnetic flux penetrating the plaquette.

Variation of Eq. (18) with respect to $`A_{𝐧\alpha }`$ gives the equation

$$I_{𝐧\alpha }=\frac{2e}{\mathrm{}}\underset{𝐤}{}\stackrel{~}{}\times M_{\mathrm{𝐧𝐤}}^1(\times A)_𝐤,$$
(25)
which relates the value of the superconducting current through a junction

$$I_{𝐧\alpha }=I_c\mathrm{sin}(_\alpha \phi _𝐧A_{𝐧\alpha })$$
(26)
with the vector potential of the magnetic field induced by the presence of the currents in the array. Here \[like in Eq. (24)\] $`\stackrel{~}{}\times `$ stands for $`_\beta ϵ_{\beta \alpha }\stackrel{~}{}_\beta `$, whereas $`\stackrel{~}{}_\beta `$ designates the lattice difference, analogous to the one defined by Eq. (21), but shifted in the negative direction:

$$\stackrel{~}{}_\beta X_𝐧X_𝐧X_{𝐧𝐞_\beta }.$$
(27)
On the other hand variation of Eq. (18) with respect to $`\phi _𝐧`$ gives the current conservation equation

$$(\stackrel{~}{}I)_𝐧=0,$$
(28)
where

$$(\stackrel{~}{}I)_𝐧\underset{\alpha }{}[I_{𝐧\alpha }I_{(𝐧𝐞_\alpha )\alpha }]$$
(29)
is the lattice equivalent of divergence. Eq. (28) can be alternatively obtained by application of the operator $`\stackrel{~}{}_\alpha `$ to Eq. (25). Therefore Eq. (25) and Eq. (28) \[both obtained by variation of Eq. (18)\] are not independent of each other.

The vector potential of the magnetic field created by the currents flowing in the array can be chosen purely transverse ($`\text{div}𝐀=0`$), which in terms of $`A_{𝐧\alpha }`$ corresponds to

$$\stackrel{~}{}_\alpha A_{𝐧\alpha }=0.$$
(30)
In that case Eq. (25) is reduced to

$$I_{𝐧\alpha }=\frac{2e}{\mathrm{}}\underset{𝐤}{}M_{\mathrm{𝐧𝐤}}^1(\mathrm{\Delta }_LA_\alpha )_𝐤,$$
(31)
where $`\mathrm{\Delta }_L_\beta \stackrel{~}{}_\beta _\beta `$ is the two-dimensional lattice analog of the Laplacian:

$$(\mathrm{\Delta }_LX)_𝐤=\underset{\beta }{}(X_{𝐤+𝐞_\beta }2X_𝐤+X_{𝐤𝐞_\beta }).$$
(32)
Comparison of Eq. (31) with Eq. (6) allows to find that for $`|𝐧𝐤|1`$

$$M_{\mathrm{𝐧𝐤}}^1\left(\frac{\mathrm{}}{2e}\right)^2\frac{1}{\pi \mu _0a|𝐧𝐤|},$$
(33)
whereas for $`|𝐧𝐤|1`$ the form of $`M_{\mathrm{𝐧𝐤}}^1`$ depends on the particular shape of superconducting islands .

Linearization of Eqs. (25) and their solution allows to show that when the magnetic fields of the currents in the array are taken into account, the logarithmic interaction of vortices becomes screened  at so-called magnetic field penetration length $`\mathrm{\Lambda }`$, exactly as it happens in superconducting films . When screening is relatively weak (that is when $`\mathrm{\Lambda }a`$), the value of $`\mathrm{\Lambda }`$ is given by

$$\mathrm{\Lambda }\frac{2}{\mu _0J}\left(\frac{\mathrm{}}{2e}\right)^2$$
(34)
and does not depend on the shape of superconducting islands forming the array .

Instead of considering Hamiltonian (18) as dependent on two different types of variables defined on the sites ($`\phi _𝐧`$) and on the bonds ($`A_{𝐧\alpha }`$) of the lattice, it is convenient to use a single variable, namely the gauge invariant phase difference $`\theta _{𝐧\alpha }`$ defined by Eq. (22). In terms of $`\theta _{𝐧\alpha }`$ the Hamiltonian (18) can be rewritten as

$$H=J\underset{𝐧,\alpha }{}\mathrm{cos}\theta _{𝐧\alpha }+\frac{1}{2}\underset{𝐧,𝐤}{}(\times \theta )_𝐧M_{\mathrm{𝐧𝐤}}^1(\times \theta )_𝐤,$$
(35)
variation of which with respect to $`\theta _{𝐧\alpha }`$ reproduces Eq. (25) in the form

$$I_{𝐧\alpha }=\frac{2e}{\mathrm{}}\underset{𝐤}{}\stackrel{~}{}\times M_{\mathrm{𝐧𝐤}}^1(\times \theta )_𝐤,$$
(36)
where the expression for the superconducting current

$$I_{𝐧\alpha }=I_c\mathrm{sin}\theta _{𝐧\alpha }$$
(37)
is naturally consistent with Eq. (26). As previously, the current conservation equation (28) can be obtained by application of the operator $`\stackrel{~}{}_\alpha `$ to Eq. (36).

## IV Dynamic fluctuations in array of resistively shunted junctions

The dynamic description of the same system requires to complement the Hamiltonian $`H`$ by the dissipative function $`W`$ (we assume that the array is overdamped and therefore its dynamics is purely relaxational). In the case of the array formed by SNS (superconductor - normal metal - superconductor) junctions one can describe dissipation in terms of the effective resistance shunting each junction (so-called RSJ-model). This corresponds to $`W\{\theta \}`$ of the form

$$W=\eta \underset{𝐧,\alpha }{}\left(\frac{}{t}\theta _{𝐧\alpha }\right)^2,$$
(38)
where the effective viscosity

$$\eta =\left(\frac{\mathrm{}}{2e}\right)^2\frac{1}{R}$$
(39)
is determined by the value of the shunting resistance $`R`$, which is assumed to be the same for all junctions. For $`W`$ of the form (38) the conservation of energy is achieved when the time evolution of the variables $`\theta _{𝐧\alpha }`$ is governed by the standard equations of relaxational dynamics:

$$\eta \frac{}{t}\theta _{𝐧\alpha }=\frac{H}{\theta _{𝐧\alpha }}.$$
(40)
On the other hand, Eq. (40) can be rewritten in the form (36), where the expression for the current should be replaced by

$$I_{𝐧\alpha }=I_c\mathrm{sin}\theta _{𝐧\alpha }+\frac{\mathrm{}}{2eR}\frac{}{t}\theta _{𝐧\alpha }.$$
(41)
The time derivative of $`\theta _{𝐧\alpha }`$ being proportional to the voltage, the second term in Eq. (41) can be easily identified as the normal current flowing through the junction. Consideration of purely relaxational dynamics means that we are neglecting capacitive effects and currents have to be conserved on each site of the lattice (in other words, only transverse current are allowed). This is ensured by the form of Eq. (36), substitution of which into Eq. (28) automatically leads to its fulfillment for any form of $`I_{𝐧\alpha }`$.

In presence of thermal fluctuations the right-hand side of Eq. (40) should be complemented with the random force term $`\xi _{𝐧\alpha }(t)`$:

$$\eta \frac{}{t}\theta _{𝐧\alpha }=\frac{H}{\theta _{𝐧\alpha }}+\xi _{𝐧\alpha }+f_{𝐧\alpha },$$
(42)
the correlations of which are Gaussian and satisfy

$$\xi _{𝐧\alpha }(t)\xi _{𝐤\beta }(t^{})=2\eta T\delta _{\mathrm{𝐧𝐤}}\delta _{\alpha \beta }\delta (tt^{}),$$
(43)
where $`T`$ is the temperature expressed in energy units (that is multiplied by the Boltzmann constant $`k_B`$). We also have included in the right-hand side of Eq. (42) the non-random external force $`f_{𝐧\alpha }`$ (to be discussed later).

In terms of the expression for the current the inclusion into Eq. (42) of the random term $`\xi _{𝐧\alpha }`$ corresponds to appearance in Eq. (41) of the additional (fluctuating) contribution to normal current $`\delta I_{𝐧\alpha }`$:

$$I_{𝐧\alpha }=I_c\mathrm{sin}\theta _{𝐧\alpha }+\frac{\mathrm{}}{2eR}\frac{}{t}\theta _{𝐧\alpha }+\delta I_{𝐧\alpha },$$
(44)
where $`\delta I_{𝐧\alpha }(2e/\mathrm{})\xi _{𝐧\alpha }`$. Note that since Eq. (36) describes the relation between the currents and the magnetic field induced by them, it has to remain fulfilled also when fluctuations of currents are taken into account. The validity of the current conservation equations (28) remains ensured by the form of the right-hand side of Eq. (36).

The suggestion to describe the dynamics of a resistively shunted Josephson junction array by Eqs. (44) has been put forward by Shenoy , who did not consider fluctuations of the magnetic field, that is assumed $`\theta _{𝐧\alpha }\phi _{𝐧+𝐞_\alpha }\phi _𝐧`$. In that case substitution of Eqs. (44) into current conservation equations (28) leads to dynamic equations for $`\phi _𝐧`$ with non-local effective viscosity . Quite remarkably the inclusion into consideration of the magnetic field fluctuations leads to simplification of the dynamic equations which (in terms of $`\theta _{𝐧\alpha }`$) become local. The idea that in presence of magnetic field fluctuations a resistively shunted Josephson junction array can be described by Eqs. (42), where $`\theta _{𝐧\alpha }`$ is the gauge invariant phase difference, has been introduced by Domínguez and José .

## V The fluctuation-dissipation theorem

It is well known that when the time evolution of some variables $`\{\theta \}`$ is governed by the standard Langevin equations of the form (42), the equilibrium (that is calculated for $`f_{𝐧\alpha }=0`$) correlation function

$$C_{𝐧\alpha ,𝐤\beta }(tt^{})\theta _{𝐧\alpha }(t)\theta _{𝐤\beta }(t^{})_{f=0}$$
(45)
is related with the response function

$$G_{𝐧\alpha ,𝐤\beta }(tt^{})\frac{\delta \theta _{𝐧\alpha }}{\delta f_{𝐤\beta }}|_{f=0}$$
(46)
by the fluctuation-dissipation theorem:

$$G_{𝐧\alpha ,𝐤\beta }(t)G_{𝐧\alpha ,𝐤\beta }(t)=\frac{1}{T}\frac{}{t}C_{𝐧\alpha ,𝐤\beta }(t).$$
(47)
However, in practical situation one is interested not in the response of the gauge invariant phase difference $`\theta _{𝐧\alpha }`$ to the (unspecified) conjugate force $`f_{𝐧\alpha }`$, but rather in more readily observable quantities such as conductivity, which is the response of a current to application of electric field. In situation when electric field $`𝐄(𝐫)`$ is created due to presence of ac magnetic field it is given by

$$𝐄(𝐫)=\frac{}{t}𝐀(𝐫),$$
(48)
where $`𝐀(𝐫)`$ is the vector potential defining the (total) magnetic field $`𝐁(𝐫)=\text{rot}𝐀(𝐫)`$.

In presence of the external magnetic field $`𝐁^\mathrm{e}(𝐫)=\text{rot}𝐀^e(𝐫)`$ the expression (23) describing the magnetic field energy should be replaced by the expression for the Gibbs free energy

$$F_{\mathrm{mf}}=\frac{1}{2\mu _0}d^3𝐫[𝐁^2(𝐫)2𝐁(𝐫)𝐁^\mathrm{e}(𝐫)],$$
(49)
variation of which in absence of any other terms gives $`𝐁(𝐫)=𝐁^\mathrm{e}(𝐫)`$. This leads to appearance in the Hamiltonian (35), describing the array, of the additional term

$$\mathrm{\Delta }H=\underset{𝐧,𝐤}{}(\times \theta )_𝐧M_{\mathrm{𝐧𝐤}}^1(\times A^\mathrm{e})_𝐤.$$
(50)
Here the variables

$$A_{𝐧\alpha }^\mathrm{e}=\frac{2e}{\mathrm{}}_{a𝐧}^{a(𝐧+𝐞_\alpha )}𝑑𝐫𝐀^e(𝐫),$$
(51)
defined on lattice bonds, are related to the vector potential $`𝐀^\mathrm{e}(𝐫)`$ of the external magnetic field exactly in the same way as earlier introduced variables $`A_{𝐧\alpha }`$ are related to the total vector potential $`𝐀(𝐫)`$.

The form of the correction to the Hamiltonian given by Eq. (50) corresponds to the presence in Eq. (42) of the external force term

$$f_{𝐧\alpha }=\underset{𝐤}{}\stackrel{~}{}\times M_{\mathrm{𝐧𝐤}}^1(\times A^\mathrm{e})_𝐤.$$
(52)
Comparison of Eq. (52) with Eq. (36) shows that (up to the factor of $`2e/\mathrm{}`$) $`f_{𝐧\alpha }`$ is related to $`A_{𝐧\alpha }^\mathrm{e}`$ exactly in the same way as $`I_{𝐧\alpha }`$ is related to $`\theta _{𝐧\alpha }`$. This allows to conclude that the correlation functions of the currents in the array

$$C_{𝐧\alpha ,𝐤\beta }^I(tt^{})I_{𝐧\alpha }(t)I_{𝐤\beta }(t^{})$$
(53)
and the response function

$$G_{𝐧\alpha ,𝐤\beta }^I(tt^{})\frac{\delta I_{𝐧\alpha }}{\delta (A_{𝐤\beta }^\mathrm{e})^t}|_{f=0}$$
(54)
have to satisfy the relation

$$G_{𝐧\alpha ,𝐤\beta }^I(t)G_{𝐧\alpha ,𝐤\beta }^I(t)=\frac{\mathrm{}}{2e}\frac{1}{T}\frac{}{t}C_{𝐧\alpha ,𝐤\beta }^I(t)$$
(55)
completely analogous to Eq. (47). Here $`(A_{𝐤\beta }^\mathrm{e})^t`$ is the transverse part of $`A_{𝐤\beta }^\mathrm{e}`$. As can be seen from the right-hand side of Eq. (52), the longitudinal part of $`A_{𝐤\beta }^\mathrm{e}`$ is completely decoupled from fluctuations of $`\theta _{𝐧\alpha }`$.

In terms of effective conductivity $`g_{𝐧\alpha ,𝐤\beta }^{\mathrm{eff}}`$, defined as the response function

$$g_{𝐧\alpha ,𝐤\beta }^{\mathrm{eff}}(tt^{})\frac{\delta }{\delta V_{𝐤\beta }^\mathrm{e}}I_{𝐧\alpha }|_{V^\mathrm{e}=0}$$
(56)
of the current with respect to external voltage

$$V_{𝐤\beta }^e=\frac{\mathrm{}}{2e}\frac{}{t}(A_{𝐤\beta }^\mathrm{e})^t,$$
(57)
Eq. (55) can be rewritten as

$$C_{𝐧\alpha ,𝐤\beta }^I(t)=T[g_{𝐧\alpha ,𝐤\beta }^{\mathrm{eff}}(t)+g_{𝐧\alpha ,𝐤\beta }^{\mathrm{eff}}(t)].$$
(58)
For wavevectors small in comparison with $`a^1`$, the variables $`I_{𝐧\alpha }`$ can be identified with $`a𝐢(𝐱)`$, whereas $`V_{𝐧\alpha }^\mathrm{e}`$ with $`a𝐄_{}^e(𝐱)`$, where $`𝐄_{}^\mathrm{e}(𝐱)`$ is the projection of the external electric field $`𝐄^\mathrm{e}(𝐱)`$ on the plane $`x_3=0`$ (here and below we assume that $`𝐄_{}^\mathrm{e}(𝐪,\omega )`$ is transverse). This allows to rewrite Eq. (58) as

$$𝐢^t𝐢^t_{𝐪\omega }=2T\text{Re}\left[g^{\mathrm{eff}}(𝐪,\omega )\right],$$
(59)
where the Fourier transform $`g^{\mathrm{eff}}(𝐪,\omega )`$ of the effective conductivity is the coefficient of proportionality in the relation

$$𝐢(𝐪,\omega )=g^{\mathrm{eff}}(𝐪,\omega )𝐄_{}^\mathrm{e}(𝐪,\omega ).$$
(60)
One should distinguish between $`g^{\mathrm{eff}}(𝐪,\omega )`$ and (also momentum and frequency dependent) sheet conductance $`g_{\mathrm{}}(𝐪,\omega )`$, which is defined as the coefficient of proportionality between $`𝐢(𝐪,\omega )`$ and total electric field $`𝐄_{}(𝐪,\omega )`$:

$$𝐢(𝐪,\omega )=g_{\mathrm{}}𝐪,\omega )𝐄_{}(𝐪,\omega ).$$
(61)
The form of the current induced correction to $`𝐄_{}^\mathrm{e}(𝐪,\omega )`$ can be easily found with the help of Eq. (6) and Eq. (48), which lead to

$$𝐄_{}(𝐪,\omega )=𝐄_{}^\mathrm{e}(𝐪,\omega )+i\omega \frac{\mu _0}{2q}𝐢(𝐪,\omega ).$$
(62)
Substitution of Eq. (61) into Eq. (62) then gives

$$𝐄_{}(𝐪,\omega )=\frac{1}{1i\omega \frac{\mu _0}{2q}g_{\mathrm{}}(𝐪,\omega )}𝐄_{}^\mathrm{e}(𝐪,\omega )$$
(63)
and, accordingly,

$$g^{\mathrm{eff}}(𝐪,\omega )=\frac{1}{i\omega \frac{\mu _0}{2q}+Z_{\mathrm{}}(𝐪,\omega )},$$
(64)
where $`Z_{\mathrm{}}(𝐪,\omega )1/g_{\mathrm{}}(𝐪,\omega )`$ is momentum and frequency dependent sheet impedance.

The form of Eq. (64) suggests that the response of a current in a two-dimensional system to external electric field is the same as if the proper sheet impedance of a system $`Z_{\mathrm{}}(𝐪,\omega )`$ has been connected in series with the other contribution, which can be considered as the effective impedance of the empty space surrounding this system. This additional contribution is purely inductive and corresponds to

$$L_s(q)=\frac{\mu _0}{2q}.$$
(65)
In the case of a superconducting system in the low frequency limit $`Z_{\mathrm{}}(𝐪,\omega )i\omega L_{\mathrm{}}`$, where $`L_{\mathrm{}}`$ is the effective sheet inductance, substitution of which into Eq. (63) allows to rewrite it as

$$𝐄_{}(𝐪,\omega )=\frac{\mathrm{\Lambda }q}{1+\mathrm{\Lambda }q}𝐄_{}^\mathrm{e}(𝐪,\omega )$$
(66)
where

$$\mathrm{\Lambda }=\frac{2L_{\mathrm{}}}{\mu _0}$$
(67)
is the (two-dimensional) magnetic field penetration length  already discussed in Sec. 3. Since in the considered system electric field appears only due to presence of ac magnetic field, the same length describes as well the screening of the electric field.

The form of the fluctuation-dissipation theorem obtained after substitution of Eq. (64) into Eq. (59)

$$𝐢^t𝐢^t_{𝐪\omega }=2T\text{Re}\frac{1}{i\omega L_s(q)+Z_{\mathrm{}}(q,\omega )}$$
(68)
being completely independent of the details of the structure of the particular system used for its derivation, one can expect it also to be valid for other two-dimensional superconducting (or simply conducting) systems, in particular superconducting films.

## VI Results and discussion

Substitution of Eq. (68) into Eq. (16) gives the expression for the flux noise power spectrum

$$S(\omega )=\frac{\mu _0^2T}{2}\frac{d^2𝐪}{(2\pi )^2}F^2(q)\text{Re}\left[\frac{1}{i\omega L_s(q)+Z_{\mathrm{}}(q,\omega )}\right]$$
(69)
which is the central result of this work. It allows instead of constructing special theories explaining frequency dependence of $`S(\omega )`$ in different regimes to relate it with already known properties of $`Z_{\mathrm{}}(q,\omega )`$.

It is of interest to compare Eq. (69) with the expression for the quantity \[the correction to frequency dependent mutual impedance $`Z_m(\omega )`$\], which is measured in the framework of the two coil method and defined as the ratio of the measured voltage to the driving current. In our notation this expression  can be rewritten as

$$\delta Z_m(\omega )=\omega ^2\frac{\mu _0^2}{4}\frac{d^2𝐪}{(2\pi )^2}\frac{F(q)\stackrel{~}{F}(q)}{i\omega L_s(q)+Z_{\mathrm{}}(q,\omega )},$$
(70)
where $`\stackrel{~}{F}(q)`$ is the geometrical factor analogous to $`F(q)`$, but describing the driving coil. Note that applicability of Eq. (69) and Eq. (70) for description of a Josephson junction array requires the lattice constant of the array to be much smaller then the geometrical parameters of the coil(s) ($`r`$ and $`h`$) entering the expressions for $`F(q)`$ \[and $`\stackrel{~}{F}(q)`$\].

Comparison of Eq. (70) with Eq. (69) shows that the real part of $`\delta Z_m(\omega )`$ is determined by the same (real) component of $`g^{\mathrm{eff}}(q,\omega )`$ as the noise spectrum, so in absence of a difference between $`F(q)`$ and $`\stackrel{~}{F}(q)`$ the noise spectrum $`S(\omega )`$ would be simply proportional to $`\text{Re}[\delta Z_m(\omega )]`$:

$$S(\omega )=\frac{2T}{\omega ^2}\text{Re}[\delta Z_m(\omega )].$$
(71)
In the case of the simplest coil (a single circular loop of radius $`r`$) Eq. (69) is reduced to

$`S(\omega )=\pi \mu _0^2Tr^2{\displaystyle _0^{\mathrm{}}}`$ $`dq`$ $`\{{\displaystyle \frac{J_1^2(rq)}{q}}\mathrm{exp}(2hq)`$ (72)
$`\times `$ $`\text{Re}\left[{\displaystyle \frac{1}{i\omega L_s(q)+Z_{\mathrm{}}(q,\omega )}}\right]\}`$ (73)
The analogous equation derived by Kim and Minnhagen  differs from Eq. (72) basically (i) by the absence of the term $`L_s(q)`$ and (ii) by the presence in the integrand of the additional factor $`q^4`$. The former of the two discrepancies is rather natural, since in Ref.  the screening effects have not been taken into account, whereas the latter we believe to be the consequence of the incorrect calculation of the magnetic field produced by currents in the array.

In Ref.  this magnetic field is calculated as the sum of magnetic fields produced by current loops associated with lattice plaquettes, whereas the magnitude of a current in each loop is assumed to be given the directed sum of the currents in the junctions $`I_{𝐧\alpha }`$ along the perimeter of a plaquette (in the present work this sum is denoted $`\times I`$). In such procedure the current of each junction is counted twice (as giving contributions to the loop currents associated with the two neighboring lattice plaquettes) but with opposite signs, so these two contributions almost cancel each other, which leads to appearance of the additional $`q`$-dependent factor in comparison with our Eq. (72). In a more consistent implementation of this approach the values of loop currents associated with arrays plaquettes should be chosen in such a way that the value of the current on each junction is given by the difference of the loop currents associated with the two neighboring plaquettes. These so called mesh currents  $`I_𝐧^\mathrm{m}`$ are related with the currents on the junctions $`I_{𝐧\alpha }`$ as

$$I_𝐧^\mathrm{m}\mathrm{\Delta }_L^1(\times I)_𝐧$$
(74)
which explains the appearance of the additional factor of $`q^4`$ in calculation which uses $`(\times I)_𝐧`$ instead of $`I_𝐧^\mathrm{m}`$.

The main contribution to the integral in Eq. (72) is coming from the region

$$0<q1/h,$$
(75)
so in the situation when $`h`$ exceeds all the microscopic scales responsible for the $`q`$ dependence of $`Z_{\mathrm{}}(q,\omega )`$, one can replace in Eq. (72) $`Z_{\mathrm{}}(q,\omega )`$ by

$$Z_{\mathrm{}}(\omega )i\omega L_{\mathrm{}}(\omega )+R_{\mathrm{}}(\omega )=\underset{q0}{lim}Z_{\mathrm{}}(q,\omega )$$
(76)
and use $`S(\omega )`$ to extract information about $`L_{\mathrm{}}(\omega )`$ and $`R_{\mathrm{}}(\omega )`$.

Let us first discuss the limit when the effects of screening can be neglected. This is possible if in the essential part of the interval (75) one can neglect $`L_s(q)`$ in comparison with $`L_{\mathrm{}}(\omega )`$ and therefore requires $`h\mathrm{\Lambda }`$ (and not $`h\mathrm{\Lambda }`$ as has been claimed in Ref.  and Ref. ). In that case Eq. (72) is reduced to

$$S(\omega )=\pi \mu _0^2Tr^2Y\left(h/r\right)\frac{R_{\mathrm{}}(\omega )}{\omega ^2L_{\mathrm{}}^2(\omega )+R_{\mathrm{}}^2(\omega )},$$
(77)
where

$`Y\left(h/r\right)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑q{\displaystyle \frac{J_1^2(rq)}{q}}\mathrm{exp}\left(2hq\right)`$ (78)
$``$ $`\{\begin{array}{cc}\frac{1}{2}\hfill & \text{for }h/r1\hfill \\ \frac{1}{16}\left(\frac{r}{h}\right)^2\hfill & \text{for }h/r1.\hfill \end{array}`$ (81)
This means that with increase of $`h`$ the amplitude of the noise has first to change very slowly and then to decay proportionally to $`1/h^2`$. The experimental results of Festin et al. , obtained on thin YBCO film at $`h/r1`$, are compatible with $`1/h^2`$ dependence even better than with $`1/h^3`$ dependence predicted by Kim and Minnhagen .

In the opposite limit of strong screening ($`\mathrm{\Lambda }h`$) $`L_{\mathrm{}}(\omega )`$ is negligibly small in comparison with $`L_s(\omega )`$ and $`S(\omega )`$ depends only on $`R_{\mathrm{}}(\omega )`$, but not on $`L_{\mathrm{}}(\omega )`$. In particular, for $`hr`$ integration in Eq. (72) gives

$$S(\omega )\{\begin{array}{cc}\frac{\pi \mu _0rT}{\omega }\hfill & \text{for }\omega \frac{R_{\mathrm{}}(\omega )}{\mu _0h}\hfill \\ \frac{2r}{h}\frac{R_{\mathrm{}}(\omega )T}{\omega ^2}\hfill & \text{for }\omega \frac{R_{\mathrm{}}(\omega )}{\mu _0h}.\hfill \end{array}$$
(82)
Note that for $`\mathrm{\Lambda }hr`$ and small $`\omega `$ one obtains the dependence $`S(\omega )1/\omega `$, which is observed in noise measurements in Josephson junction arrays  over four decades in frequency.

The idea that this particular dependence appears because SQIUD integrates contributions from wide interval of wavevectors has been put forward by Wagenblast and Fazio . Our analysis has shown that in order to obtain the $`1/\omega `$ dependence over four decades in frequency one should have $`r/h10^4`$, which in the experiments of Refs.  definitely has not been fulfilled. Thus the origin of the $`1/\omega `$ dependence of the flux noise power spectrum observed experimentally in Josephson junction arrays (in particular, at tempratures lower then the expected temperature of the Berezinskii-Kosterlitz-Thouless (BKT) phase transition ) still remains to be elucidated.

Substitution of the results of Ambegaokar et al.  for vortex pair contribution to impedance into Eq. (69) or its simplified analogs allows to reproduce the dependence $`S(\omega )1/\omega `$ only exactly at the temperature of the BKT transition. On the other hand, the analogous frequency dependence of the noise spectrum observed in thin YBCO films  is ascribed to vortex hopping between neighboring pinning centers, a finite vortex concentration being associated with presence of a residual magnetic field. The same mechanism may be responsible for the results of experiments  on Josephson junction arrays.

## Acknowledgments

The author is grateful to V. V. Lebedev for useful discussion. This work has been supported by the Program ”Quantum Macrophysics” of the Russian Academy of Sciences, by the Program ”Scientific Schools of the Russian Federation” (grant No. 00-15-96747), by the Swiss National Science Foundation, and by the Netherlands Organization for Scientific Research (NWO) in the framework of Russian-Dutch Cooperation Program.


format: text

text:
# Molecular Dynamics of Comminution in Ball Mills

## 1 Introduction

In the past years molecular dynamics simulations of granular systems – such as sand, fertilizer, grain and others – have been developed towards a reliable tool for an investigation of granular systems. Unlike a few years ago, when physicists started to investigate granular assemblies consisting of just a few hundred spherical granular particles in two dimensions, e.g. HW ; CS ; Herrmann ; Taguchi ; Walton ; WaltonBraun ; BuchholtzPoeschel:1997 , and less than 100 particles in three dimensions, e.g. GHPS , today, owing to the rapid evolution of computer facilities, we are able to simulate complex systems of many thousands of particles in two and three dimensions even accounting for non-spherical grain features.

Due to this progress molecular dynamics of granular systems by now can be applied to the simulation of technologically important processes striving for a better understanding of relevant details of the process and, henceforth, for an optimization of industrial technologies (e.g. Kohring:1993 ; SEP ; SEP1 ; Kohring ; YokoyamaTamuraJimbo:1994 ; YokoyamaTamuraUsuiJimbo:1996 ; SongfackAgrawalaMishraRajamani:1993 ). Molecular simulation techniques, therefore, might develop to a prospective engineering tool partly replacing costly laboratory experiments by “computer experiments”.

In the present paper we apply the method of molecular dynamics to the simulation of autogenous dry comminution in a tumbling ball mill ballmill . To this end we devised a novel algorithm which accounts for the fragmentation of particles. First we will demonstrate that experimentally known results, in particular those by Rothkegel and Rolf Rothkegel:1992 ; Rolf , can be reproduced up to a good degree of accuracy. We will show that the efficiency of comminution in ball mills is mainly determined by the presence of force chains. Our result will hint at how to improve the efficiency of a milling machinery widely applied in industry.

## 2 The Simulation Model

### 2.1 Molecular Dynamics

For the simulation we assumed pairwise interaction forces and applied the model by Cundall and Strack CS  and Haff and WernerHW : A pair of two-dimensional spherical particles $`i`$ and $`j`$ interacts through the force

$$\stackrel{}{F}_{ij}=F_{ij}^N\stackrel{}{n}+F_{ij}^T\stackrel{}{t},$$
(1)
if their distance is smaller than the sum of their radii $`R_i+R_j\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|=\xi _{ij}>0`$. The unit vectors in normal and tangential direction are defined as

$`\stackrel{}{n}`$ $`=`$ $`\stackrel{}{r}_{ij}/\left|\stackrel{}{r}_{ij}\right|`$ (2)
$`\stackrel{}{t}`$ $`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{01}{10}}\right)\stackrel{}{r}_{ij}/\left|\stackrel{}{r}_{ij}\right|,`$ (3)
with $`\stackrel{}{r}_{ij}\stackrel{}{r}_i\stackrel{}{r}_j`$. The forces in normal and tangential direction are respectively given by

$`F_{ij}^N`$ $`=`$ $`Y\xi _{ij}m_{ij}^{\mathrm{eff}}\gamma _N\dot{\stackrel{}{r}}_{ij}\stackrel{}{n}`$ (4)
$`F_{ij}^T`$ $`=`$ $`\text{sign}(v_{ij}^{\mathrm{rel}})\mathrm{min}(m_{ij}^{\mathrm{eff}}\gamma _T\left|v_{ij}^{\mathrm{rel}}\right|,\mu \left|F_{ij}^N\right|)`$ (5)
with

$`v_{ij}^{\mathrm{rel}}`$ $`=`$ $`\dot{\stackrel{}{r}}_{ij}\stackrel{}{t}+R_i\mathrm{\Omega }_i+R_j\mathrm{\Omega }_j`$ (6)
$`m_{ij}^{\mathrm{eff}}`$ $`=`$ $`{\displaystyle \frac{m_im_j}{m_i+m_j}}.`$ (7)
Here $`R_i`$, $`m_i`$, $`\stackrel{}{r}_i`$, $`\dot{\stackrel{}{r}}_i`$ and $`\mathrm{\Omega }_i`$ are respectively the radius, the mass, the position, the velocity and the rotation velocity of the $`i`$th particle. $`Y=810^6`$ g s-2 is the elastic constant and $`\gamma _N=800`$ s-1, $`\gamma _T=3000`$ s-1 are the parameters of damping in normal and tangential direction. The Coulomb friction constant was assumed to be $`\mu =0.5`$. These empirical values which we used throughout all simulations have proven to render realistic behaviour for a typical granular material.

The normal force (4) is composed of an elastic repulsive part and a dissipative part which acts against the direction of motion. For a collision of two-dimensional spheres, i.e. disks, the Hertz contact law $`F_{ij}^N\xi _{ij}^{3/2}`$ Hertz:1882 ; Landau:1965  reduces to Engel

$$\xi _{ij}F_{ij}^N\left(\frac{2}{3}+\mathrm{ln}\frac{4(R_i+R_j)E_r}{F_{ij}^N}\right),$$
(8)
where $`E_r`$ is the reduced elastic module, i.e. a material constant. Equation (8) provides a relation $`\xi =\xi \left(F^N\right)`$, however, in the molecular dynamics simulations we need the inverse $`F^N=F^N(\xi )`$. To calculate the forces due to Eq. (8) one would have to invert (8) numerically for each particle contact in each time step, i.e., to solve a transcendental equation. Alternatively one could tabulate the function $`F^N(\xi )`$ but due to large force gradients, in particular in the instant of contact, we have to simulate with double precision accuracy. Therefore, the table would be extremely large. Apart from the small logarithmic term $`F_{ij}^N\mathrm{ln}F_{ij}^N`$ in Eq. (8) one finds that the force is proportional to the compression $`\xi _{ij}`$. Therefore, in our simulations we used the linear law (4) which is a good approximation in the force interval of interest.

Equation (6) describes the relative velocity of the particle surfaces at the point of contact which results from both the tangential part of the relative particle velocity and the velocity of particle spin.

The Coulomb friction law is taken into account in Eq. (5): If the tangential force of two colliding particles exceeds $`\mu `$ times the normal force the particles slide upon each other feeling constant friction. In this way the Coulomb law formulates a maximum transferable shear force.

As numerical integration scheme we used the Gear predictor-corrector method of sixth order (e.g. AllenTildesley ).

Let us remark that there exist various models with diverse descriptions of the interaction forces, e.g. WolfKraefte . In three-dimensional calculations it is crucial that the dissipative term in (4) is proportional to $`\sqrt{\xi }\dot{\xi }`$ BSHP  (in viscoelastic approximation). In three dimensions the used term $`\dot{\xi }`$ in Eq. (4) yields wrong results WolfKraefte ; rospap .

Another class of molecular dynamics are event-driven simulations, e.g. Rapaport  which require much less numerical effort than the method described above (force-driven Molecular Dynamics). In event-driven simulations one does not integrate the equations of motion explicitly but rather expresses the loss of mechanical energy in normal and tangential direction due to a collision via coefficients of restitution $`ϵ^N`$ and $`ϵ^T`$ both of which can be determined from the material properties of the particles BSHP ; SchwagerTP2 . Besides some conceptional difficulties which can be overcome by applying numerical tricks luding  these algorithms are only applicable in case the duration of collisions is negligible in comparison with the mean free flight time, i.e. the time a particle moves without interaction. In this limit collisions can be assumed as instantaneous events. Moreover, this implies that three-particle interactions are very rare events. However, in the system under study, namely the ball mill, this premise does not hold because the majority of particles permanently is in close contact with neighbouring particles.

### 2.2 Fragmentation Probability of Grains

It can be observed in experiments that different particles of identical size and material vary with respect to their tensile strengths. That means when compressing particles of a sample by applying a fixed force only a fraction of the sample actually will break. For a sufficiently large sample this fraction can be rephrased in terms of a fracture probability.

The reason for this different behavior can be explained by the existence of flaws which provide sites for stress concentration and the initiation of a crack which, subsequently, propagates through the material. These flaws are distributed throughout the whole volume, however, several investigations have revealed that surface flaws activated by high tensile stress play the dominant role for the initiation of a crack Weichert:1992 . Hence, the fracture probability is intimately related to the statistical distribution of surface flaws and the resulting stress distribution over the particle’s surface.

Starting from some plausible assumptions concerning the flaw distribution and the behaviour of the material the fracture probability $`P`$ of a single particle can be derived employing statistical reasoning. A particle of radius $`R`$, subjected to a force which stores the specific elastic energy $`W_m=W/m`$, will break with probability Weichert:1988

$$P(R,W_m)1\mathrm{exp}\left(cR^2W_m^z\right).$$
(9)
Here, $`c`$ and $`z`$ are material constants. The original derivation (considering rods of a brittle material) goes back to Weibull Weibull:1939 . This simple law has proven in practice by fitting various experimental data. The constants $`c`$ and $`z`$ can be extracted from measurements when plotting $`[\mathrm{ln}\mathrm{ln}(1P)^1]`$ vs. $`[\mathrm{ln}W_m]`$. Typical values of $`z`$ are found to lie in the range $`1.5,\mathrm{},2.5`$.

Notice that in agreement with experimental results the term $`R^2`$ in the exponent (for fixed $`W_m`$) in general predicts an anticorrelation between particle size and resistance to breakage. Clearly, this fact is rooted in the diminished probability to find a sufficiently large flaw on a smaller surface.

For spherical particles (or, more generally, for particles having convex surface), assuming linear elastic behaviour, Hertz-theory Hertz:1882 ; Landau:1965  can be employed to derive a relation between the elastic energy $`W`$ stored in a sphere and the exerted repulsive contact force $`F`$

$$W=\frac{2}{5}hF=\frac{2}{5}\left(\frac{D^2F^5}{R}\right)^{\frac{1}{3}},$$
(10)
where $`h`$ quantifies the particle’s deformation and $`D`$ is another material constant (expressed by the Poisson ratio and constant of elasticity) Hertz:1882 ; Landau:1965 . It should be kept in mind that formula (10) remains valid in the linear elastic regime only.

Applied to our simulation of ball mills we make use of these laws in the following way: Consider the situation sketched in Fig. 1: A particle placed in the centre is compressed by four impacting particles.

The highest compressive stress occurs around the contact point; a fact well known to engineers BernotatSchoenert:1988  and also found in numerical simulations KunHerrmann:1996 ; KunHerrmann:1996a . Hence, the fracture probability of the centre particle is linked to the probability of finding a flaw of sufficient size in the very vicinity of any of the four contact regions. Thus, the action of different impacting particles can be regarded as independent stress sources. That means we consider the elastic energy stored by each impacting particle separately, calculate the related fracture probability according to formulae (9) and (10) as a function of the particle’s size and the maximum of all acting normal forces $`F^N`$

$$P(R,F^N)1\mathrm{exp}\left(\stackrel{~}{c}R^{\left(2z\frac{10}{3}\right)}\left(F^N\right)^{\frac{5z}{3}}\right),$$
(11)
and actually decide the question of breakage by subsequently drawing a random number for each of the contact sites.

The constant $`\stackrel{~}{c}`$ is related to the empirical constant $`c`$ in Eq. (9); the precise connection is elaborated by inserting Eq. (10) into (9) and expressing $`m`$ by the particle radius $`mR^3`$.

Up to this point we discussed the fracture probability for three-dimensional spheres. In our simulations we modelled two-dimensional spheres, i.e. circular discs of unique thickness $`L`$. The fracture probability corresponding to Eq. (11) can be obtained in analogy to (9) yielding

$$P(R,W_m)1\mathrm{exp}\left(cRLW_m^z\right).$$
(12)
Note that the surface of a sphere $`R^2`$ has to be replaced by the surface of a cylinder $`RL`$ (excluding top and bottom wall). Correspondingly Eq. (10) turns into $`h=YF^N`$ and, hence, the elastically stored energy for spheres has to be replaced by

$$W=\frac{1}{2}hF^N=\frac{(F^N)^2}{2Y}(F^N)^2.$$
(13)
Taking into account that the mass of a cylinder is proportional $`R^2`$, one obtains

$$P(R,F^N)1\mathrm{exp}\left(\stackrel{~}{c}R\left[\frac{(F^N)^2}{R^2}\right]^z\right).$$
(14)
In section 3 we will motivate the choice $`z=2`$. With this specification one finally arrives at the following expression for the fracture probability

$$P(R,F^N)1\mathrm{exp}\left(\stackrel{~}{c}R^3\left(F^N\right)^4\right).$$
(15)
## 3 Fragment size distribution

Once a flaw has been activated it forms an initial crack which rapidly propagates through the particle. Typical crack velocities in the range of $`1500m/s`$ have been measured Weichert:1992 . The formation of a variety of differently sized fragments occurs by virtue of branching cascades. Branching is governed by a balance between the energy release rate $`G`$ and the so-called crack resistivity $`B`$ Weichert:1992  the latter accounting for the creation of new surfaces. Because of the extremely short duration of the rupture process one can safely neglect external energetic contributions to this energy balance. The energy release rate $`G`$ depends on the crack velocity and the crack length whereas the crack resistivity $`B`$ depends on the crack velocity and on the number of branches. Due to a maximal crack velocity, i.e. the speed of sound, the balance between both terms requires the formation of new cracks. From these considerations it can be shown Weichert:1992  that the fragment size distribution (cumulative mass distribution) $`Q`$ is a function of the product $`R_fW_m`$, where $`R_f`$ denotes the fragment size, i.e. $`Q=Q(R_fW_m)`$.

Experimental evidence Weichert:1988  for a scaling law $`QR_f^\beta `$ is in agreement with the well-known empirical Schuhmann law Schuhmann:1940

$$Q\left(\frac{R_f}{k}\right)^\beta ,$$
(16)
which itself can be regarded as an approximation of the Rosin-Rammler law RosinRammler:1933

$$Q1\mathrm{exp}\left[\left(\frac{R_f}{k}\right)^\beta \right].$$
(17)
The variable $`k`$ has the dimension of a length and, thus, can only be identified with the size of the original particle.

A rigorous derivation of Eq. (17) – compatible with an exponent $`\beta =1`$ – starting from rather mild assumption and applying Poisson statistics was performed by Gilvary Gilvary:1961  (see also GilvaryBergstrom:1961 ).

It has to be mentioned that the derived laws Eqs. (16,17) describe the distribution of fragments only below a certain size (endoclastic vs. exoclastic distribution Gilvary:1961 ). For the largest fragments an equivalently well–defined quantitative statement does not exist. The deviation between theoretical predictions and experimental observations typically occurs for grains which together collect about 75–90% of the total mass Weichert:1992 . Moreover, some diverging theoretical mean values can be understood from the fact that the assumption of equidistributed flaws on all length scales neglects the effect of flaw depletion applying to tiny fragments. Flaw depletion thus sets a lower bound to the range of applicability. Insofar the size of fragments which is satisfactorily described by the distributions Eq. (16) and Eq. (17) is restricted to the intermediate range.

When fitting both the more refined Rosin–Rammler law Eq. (17) and the simplified Schumann law Eq. (16) to diverse experimental data it was found that in some situations the Rosin-Rammler law proved superior Weichert:1992  (mill products GilvaryBergstrom:1961 ) whereas in other cases (specimen fractured in gelatin GilvaryBergstrom:1961 ) the Schuhmann law was even more adequate.

Finally, we mention the population balance model Herbstetal:1981  which assumes that the fragment size distribution can be normalized to the initial particle size, i.e. $`Q=Q(R_f/R)`$. This assumption is based on empirical observations. Moreover, its central assertion can also be derived when starting from a Weibull exponent $`z=2`$ and applying statistical reasoning similar to the one sketched above.

In our simulations we employed the Rosin-Rammler law Eq. (17) with $`k=R`$ and an exponent $`\beta =1`$. Since experimental data have shown the mean dimension of the largest fragment to be roughly 75% of the original particle size GilvaryBergstrom:1961 , i.e. $`R_f^{max}=\frac{3}{4}R`$, the normalized fragment size distribution simply reads

$$Q\left(\frac{R_f}{R}\right)=\frac{1\mathrm{exp}\left({\displaystyle \frac{R_f}{R}}\right)}{1\mathrm{exp}\left({\displaystyle \frac{3}{4}}\right)}.$$
(18)
### 3.1 Molecular Dynamics Modelling of Fragmentation

In the present study we exclusively utilize spherical particles (disks), which means that also the fragments of a disrupted particle must be spherical. The single fracture statistics, as discussed in section 3, is put into practice by dismembering a cracking disk into two fragments. Whereas the size of one fragment is chosen in accordance with the Rosin-Rammler law the size of the second one is determined by mass conservation, i.e. for two-dimensional particles by area conservation.

In the majority of situations a breaking particle is closely surrounded by neighbouring ones. Hence, in general, it will be impossible to place the two fragments avoiding an overlap. However, since a considerable overlap corresponds to a situation of extreme compression an unrealistically strong repulsive normal force (4) together with an exceptional deformation energy would be the consequence. We bypassed this problem by modifying the interaction between the two fragments in the following way: Up to the moment of first complete separation of two such fragments $`i`$ and $`j`$ the normal force (4) at the $`k+1`$-rst time step is replaced by

FijN={
-Yξij(k)mijeffγN(-˙ri˙rj)n if (r˙ir˙j)n<0ϵ otherwiseF_{ij}^{N}=\left\{\begin{tabular}[]{l}$Y\xi_{ij}^{*}(k)-m_{ij}^{{\rm eff}}\,\gamma_{N}\,(\dot{\vec{r}}_{i}-\dot{\vec{r}}_{j})\cdot\vec{n}$\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}if~{}~{}~{}
$(\dot{\vec{r}}_{i}-\dot{\vec{r}}_{j})\cdot\vec{n}<0$\\
$\epsilon$ ~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}otherwise\\
\end{tabular}\right. (19)
with

$$\xi _{ij}^{}(k)=\left|\stackrel{}{r}_i(k1)\stackrel{}{r}_j(k1)\right|\left|\stackrel{}{r}_i(k)\stackrel{}{r}_j(k)\right|.$$
(20)
This means that the two fragments, on the one hand, resist to further compression like standard particles and, on the other hand, experience a small repelling force $`ϵ`$ which drives them apart. The repelling force is chosen sufficiently small in order to suppress unphysical energy gain of the system; this means that the unavoidable energy, which due to this force is pumped into the system anyhow, must be negligible in comparison with the mean kinetic energy of the grains. After the moment of first complete separation of the two fragments they also interact in the standard way, as dictated by Eq. (1). Of course, the interaction of two such fragments with all other particles is never modified. By trimming of $`ϵ`$ we could not only suppress the elastic energy gain but also evidenced that the period requiring a modified interaction force was comparatively short. The procedure of fragmentation is sketched in Fig. 2. We want to mention that multiple fragmentation, i.e. further fragmentation of the fragments is included in the model and is frequently observed in simulations.

Let us remark that our algorithm was flexible enough to include multiple fracture events: By this we mean the observation that a fragment was apt to break itself even before its first complete separation from its original twin fragment.

An animated sequence of fragmentation can be found in the Internet under http://summa.physik.hu-berlin.de/ $``$kies/mill/bm.html

Finally, we mention that the proposed algorithm is not the only way to model particle fragmentation. An alternative algorithm has been put forward by Åström and Herrmann AstroemHerrmann:1998 .

### 3.2 Specification of the Model System

As depicted in Fig. 3 the cylindrical wall of our system has been modelled by an ensemble of spheres possessing the same material constants as the grains inside the container. However, the motion of wall spheres is not affected by impacts but, instead, strictly governed by the continuous rotation of the cylinder. Moreover, in spite of identical elastic and dissipative constants the wall particles were considered to be much harder so they never would break.

A suitable size distribution of wall particles generates the effect of surface roughness. This type of modelling has been used in several numerical simulations and is also used in experiments to guarantee well-defined boundary properties. The wall particles have been placed along the circumference of the cylinder. Just like in realistic ball mills we evenly distributed 16 toolbars along the inner periphery of the cylinder which serve to lift particles and thus prevent them from sliding downhill. We composed the toolbars again from a collection of rigid spheres. As reported in literature BernotatSchoenert:1988  the maximum milling efficiency is achieved for a filling degree, i.e. the ratio of material and container volume, of ca. 40%.

The simulation results to be presented in the next section have been done with a system consiting in average of about 1000 particles. Due to fragmentation the number of particles was not constant, i.e., the particle number fluctuates. The container was modeled by about 500 wall particles as described above.

Initially the spheres which model the grist were randomly positioned inside the container avoiding any overlap. This artificial initial condition has to be relaxed by evolution of the system. In this way we achieved a realistic configuration for any chosen rotation velocity. Obviously the typical asymptotic configuration depends on the rotation frequency $`\nu `$, therefore, we started our numerical monitoring only after transients had died out. In Fig. 4 we present snapshots of typical configurations depending on four different rotation frequencies.

In technologically relevant ball mills there is a continuous transport of raw material into the mill and of milled material out of the device. In most industrial applications this exchange of raw material and final product is accomplished by axial transport (z-direction) of material. Of course, with our two-dimensional model we cannot account for axial transport, therefore, we simulated the material exchange in continuous operation mode in the following way: Whenever fragments were generated whose size dropped below a minimum we removed them from the container. After exact bookkeeping of the mass of removed “dust” particles we reinserted a big particle when the accumulated dust mass exceeded an upper threshold, of course, obeying mass conservation. In this way we conserved the average filling level throughout the simulation and, thus, could reach a steady state of the continuous operation mode.

## 4 Simulation Results

### 4.1 Mass distributions

As a first result we studied the frequency dependent fragment size distribution. This distribution is not necessarily equal to the distribution of section 3, which was the result of single fracture experiments.

It is known Kolmogorov:1941 ; Halmos:1944 ; Epstein:1947  that for an arbitrary fracture mechanism, which hierarchically is applied on all length scales, one will always generate a log-normal distribution of fragment sizes; this assertion is a pure consequence of the central limit theorem. However, in our modelling the ball mill is operated as an open system, i.e. there is a particle flow into and out of the mill. Therefore, analyzing the asymptotic size distribution is also of theoretical interest.

Furthermore, the explicit knowledge of the distribution plays also a practical role, because from that knowledge one can decide whether a multi-level fragmentation process is to be preferred to a single-level one or vice versa. A multi-level process can be either a distribution of the comminution process to several ball mills, where each is filled with granular material of a specified narrow size distribution, or a spatial separation of material of different size in one mill.

Fig. 5 shows the differential fragment size distribution, i.e. the mass fraction of all particles with size within a certain interval, for four different rotational frequencies. The black circles in each column correspond to measurements at different times and thus reflect the temporal variance. The time period of each simulation was 10 seconds real time (beyond transients).

Obviously the lowest curve ($`2.5`$ Hz) changes only very weakly during the simulation, which can be interpreted as a low rate of fragmentation. This curve – and with restrictions also the top curve ($`0.5`$ Hz) – is determined by the initialization. Unfortunately, the middle curves are not significantly different, which mainly is a consequence of the short simulation time (10 seconds). To answer the interesting questions (theoretical distribution, multi-level process) it will be necessary to simulate considerably longer time periods.

### 4.2 Optimization of the efficiency

In the context of industrial applications the notion of efficiency mainly focuses on the aspects of maximum comminution rate and, perhaps, of power consumption. There are many free parameters which can be varied to maximize the efficiency: speed of rotation, filling degree, size distribution, cylinder diameter, grist size distribution, various mill types, etc.

In our analysis we investigated the dependence of the comminution rate on the speed of rotation. In order to compare our numerical results with experimental data we normalized quantities in the following way: The speed of rotation is normalized to the velocity $`n_c=\sqrt{g/M}(2\pi )^1`$ for which the centrifugal force balances the gravitational force whereas the comminution rate is normalized to its measured maximum. The result of our simulation is plotted in Fig. 6 using filled circles and dashed line. Clearly, a maximum occurs around $`nn_c`$ whereas the rate rapidly diminishes for larger or smaller velocities.

The full line in Fig. 6 depicts a similar curve following from experiments. The curve shows the relative power consumption as a function of the normalized speed of rotation. Qualitatively both curves are in good agreement; systematic deviations can be explained by the following observation: Strictly speaking both curves relate to different physical quantities. The fraction rate only accounts for the energy dissipated through fractures while the experimental curve includes all dissipative processes, for instance friction between grains, frictions between the grains and the wall, dissipative impacts, heating, etc.). As a consequence, the curve for the simulation should be found below the experimental curve. Moreover, from the same reasoning it is also clear that observed deviations are minor if the main portion of the overall dissipated energy results from the fractures, i.e. around the maximum.

### 4.3 Preferred fragmentation locations

The series of pictures in Fig. 7 illustrates the spatial distribution of fragmentation locations for four different rotation velocities. We restricted our presentation to the most important lower right segment of the ball mill (clockwise rotation). To compute this distribution the coordinates of breaking particles were stored during the simulation. After the simulation the fraction of all fragmentation events occurring at places inside little cells (coarse graining) was computed. In Fig. 7 grey values code the frequency: dark means high frequency and light low frequency respectively.

In the bottom right figure one finds a few fragmentation events near the cylinder wall. The restriction to medium-grey values results from the fact, that the rare events are equidistributed across the inner perimeter.

The predominance of rather light grey (or even white) pixels in the rest of the figures hints at regions where fragmentation events are rare as compared to the sparse medium-grey to dark pixels which indicate regions of high fracture probability. The observation that most fragmentation events occur deep inside the material and not near the surface, which is heavily agitated by impacts, is, at first sight, somewhat counter-intuitive. An explanation for this fact will follow from a closer analysis of the local pressure concentration in the next section.

## 5 Spatial distribution of pressure in ball mills

In this chapter we will dwell on the physical origin of the unusual pressure distribution pressure  measured in experiments.

As can be concluded from Figs. 7 and 8, the majority of fragmentation events does not occur near the free surface of the granular material but deep inside the material near to the walls. This strange behavior was observed already in experiments by Rothkegel and Rolf Rothkegel:1992 ; Rolf  in two-dimensional model mills. They used instrumented balls which emitted a flash of light whenever the strain exceeded a fixed threshold. The positions of recorded flashes were digitized and from the statistics of these data they inferred the stationary strain distribution inside the ball mill. For various threshold values they always found the peak strain deep inside the material (see Fig. 8). This was surprising because inside the material relative grain velocities are rather low in comparison with relative velocities at the surface.

This experimental finding, together with its validation by our simulation (see the previous section), calls for an explanation. We will solve this challenging, so far open question using the powerful possibilities of molecular dynamics, namely the option to monitor physical quantities which are hardly accessible through experiments.

In the following simulations we ignored fragmentation events – which are not relevant for our explanation – simply for the sake of computational efficiency, so we could afford longer simulation time. This time the radius of the two-dimensional ball mill was set to $`M=4`$ cm. The mill was filled with $`N=800`$ spherical grains with radii equally distributed in the interval $`[0,05,0,11]`$ cm. We performed simulations with three different rotation velocities $`\mathrm{\Omega }_\mathrm{I}=2`$ Hz, $`\mathrm{\Omega }_{\mathrm{II}}=10`$ Hz and $`\mathrm{\Omega }_{\mathrm{III}}=19`$ Hz.

Fig. 9 shows snapshots of our performed simulations. The pictures in the top row are taken form a simulation with rotation velocity $`\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{I}`$. The velocity is sufficiently enough to observe a continuous flow at the surface. The surface shape is similar to a plane. The regime with discontinuous flow (for lower rotation velocities) will not be investigated in this paper, because it is not relevant for technical purposes. For a detailed discussion of the interesting observations at the transition between continuous and stick-slip flow we refer to the work by Rajchenbach Rajchenbach:1990 .

The picture in the middle row of Fig. 9 depicts snapshots of simulation with $`\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{II}}`$. The surface of the material no longer shapes a plain but the grains fall down a steep slope hitting a flat surface.

The bottom row of Fig. 9 exhibits the snapshot of a simulation with the highest rotation velocity $`\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{III}}`$. Here the grains are carried away with the fast rotation of the cylinder. Beyond a certain point they are likely to loose contact with the wall and follow a free parabola before impacting back down on the surface.

Moreover, the snapshots of Fig. 9 provide an intuitive understanding, what mechanism is responsible for the pressure distribution mentioned above. Again, grey values in these pictures codes the local pressure $`P_i`$ acting on grain $`i`$ given by

$$P_i=\underset{j}{}F_{ij}^N.$$
(21)
The index $`j`$ is running over all particles which are in contact with grain $`i`$. Dark circles denote high pressure, light spheres low pressure.

Obviously many of the heavier stressed particles are deep inside the material, which is in good agreement with the experimental observation by Rothkegel and Rolf Rothkegel:1992 ; Rolf . Most importantly, the pictures reveal that grains experiencing peak stress often form linear structures. In the following we focus on the development and the properties of such force chains. By numerical analysis we will prove that those force chains are the crucial physical reason for the observed pressure distribution.

We defined force chains by a set of three self-evident conditions: Grains $`i`$, $`j`$ and $`k`$ are considered to be members of the same force chain if:

1. Particles $`i`$, $`j`$ and $`j`$, $`k`$ are next neighbours.
2. The pressure acting on each of the grains exceeds a certain threshold ($`P_i>1000`$ g cm s-2).
3. The connecting lines between $`i`$, $`j`$ and $`j`$, $`k`$ form an angle larger than $`150^o`$, i.e. the centres of three grains almost fall on a line.

These three conditions are evaluated by a computer algorithm. In Fig. 9 all grains belonging to a force chain are graphically connected by lines.

We see that most of the harder stressed grains could be assigned to a force chain. From this we conclude that the main part of the static and dynamic pressure propagates along the force chains. The pressure acting on a highly stressed grain inside a force chain is up to 100 times larger as compared to the average pressure of the neighbouring grains not belonging to a force chain, i.e. the force distribution is strongly inhomogeneous.

The development of such force chains is not restricted to grains in a ball mill but has been observed in different granular systems by various experimentalists DuranMazoziLudingClementRajchenbach:1996 ; Behringer:1993 ; LiuNagelSchecterCoppersmithMajudarNarayanWitten:1995 ; Dantu:1967  and numerically Schoellmann:1999 ; EsipovPoeschel:1997 ; Wolf:1996 . It could be shown, that one of the granular phases is characterized by the observation of force chains EsipovPoeschel:1997 . From this one might conclude that the occurrence of force chains is an inherent feature of all granular matter, which may substantially contribute to its specific characteristics.

Fig. 10 comprises six plots which come in three pairs; the lower plot of each pair always depicts the frequency distribution of force chains whereas the related upper one presents the (statistically evaluated) maximal pressure both as a function of the length $`L`$ of a force chain (which is nothing but the number of particles belonging to the chain). These pairs are shown for three different rotation velocities varying from $`\mathrm{\Omega }_\mathrm{I}`$ (top) to $`\mathrm{\Omega }_{\mathrm{III}}`$ (bottom).

The frequency of a force chain drops exponentially with increasing length $`L`$. For the highest rotation velocity ($`\mathrm{\Omega }_{\mathrm{III}}`$) the statistics becomes unreliable for chains longer than $`L=13`$. This behaviour, which is in contrast to observations for slower measured rotations, is due to the higher energy input, which causes a decompactification of the material.

Concerning the earlier mentioned spatial strain distribution the maximal pressure measured in each force chain is of paramount interest rather than the average pressure in a force chain (see section 2.2). For lower rotation velocities $`\mathrm{\Omega }_\mathrm{I}`$ and $`\mathrm{\Omega }_{\mathrm{II}}`$ a monotonous increase of the maximal pressure with the length $`L`$ can be observed. For a discussion of these curves one has to take into account that the pairs of plots in Fig. 10 do not constitute independent quantities. Choosing a sample of $`L`$ particles from an ensemble of particles on which forces with a given distribution are acting, the maximal force in this sample will increase with increasing $`L`$. Therefore, in order to prove significance, the result in Fig. 10 has to be weighted against this purely statistical effect.

The probability to measure a maximal pressure $`P^{\mathrm{max}}`$ in the interval $`P^{\mathrm{max}}[x;x+dx]`$ in a sample of size $`L`$ is equal to the probability that all particles feel a pressure $`P<x+dx`$ minus the probability that all $`L`$ particles feel a pressure $`P<x`$.

$$p_L(P^{\mathrm{max}}[x;x+dx])=\left[p(P<x+dx)\right]^L\left[p(P<x)\right]^L$$
(22)
Assuming equally distributed pressures $`P[0:P_m]`$, $`p(P<x)`$ means

$$p(P<x)=\underset{0}{\overset{x}{}}\frac{1}{P_m}𝑑x^{}=\frac{x}{P_m}.$$
(23)
By inserting Eq. (23) into Eq. (22), one finds in first order of $`dx`$

$$p_L(P^{\mathrm{max}}[x;x+dx])=\frac{1}{\left(P_m\right)^L}Lx^{L1}dx.$$
(24)
Therefore, the average maximal pressure as a function of sample size $`L`$ is given by

$`P^{\mathrm{max}}(L)`$ $`=`$ $`{\displaystyle \frac{1}{\left(P_m\right)^L}}{\displaystyle \underset{0}{\overset{P_m}{}}}xLx^{L1}𝑑x`$ (25)
$`=`$ $`{\displaystyle \frac{1}{\left(P_m\right)^L}}{\displaystyle \frac{L}{L+1}}P_m^{L+1}`$
$`=`$ $`{\displaystyle \frac{L}{L+1}}P_m.`$
The dashed line in Fig. 10 shows this statistical effect for $`P^{\mathrm{max}}`$ with $`P_m=4800g`$ cm s-2. Obviously, these curves do not reach up to our numerical data of the top and center figure; especially for larger $`L`$ the dashed line falls far below. The numerical data points in the bottom figure (high rotation velocity) to some extent coincide with the dashed curve, which only represents the statistical effect discussed above. Hence, we can conclude that for lower velocities the maximum stress significantly concentrates within long force chains whereas for high rotation velocities a significant effect of force chains is not evident.

From experimental and numerical results (for instance LiuNagelSchecterCoppersmithMajudarNarayanWitten:1995 ; Howell:1997 ; Radjai:1996 ) it is known, that the observation of force chains is connected to a strongly inhomogeneous pressure distribution. The distribution is essentially an exponentially decreasing function with increasing pressure. A simplified expression for this distribution is given by

$$p(P)=a\mathrm{exp}\{aP\}$$
(26)
In analogy to Eq. (23) the probability to measure a pressure $`P<x`$ is given by

$$p(P<x)=\underset{0}{\overset{x}{}}a\mathrm{exp}\{ax^{}\}𝑑x^{}=1\mathrm{exp}\{ax\}.$$
(27)
Inserting this into Eq. (22) we obtain in the first order of $`dx`$

$`p_L(P^{\mathrm{max}}[x;x+dx])=`$
$`=\left(1e^{a(x+dx)}\right)^L\left(1e^{ax}\right)^L`$
$`=L(1e^{ax})^{L1}e^{ax}adx.`$ (28)
The average of the maximal pressure can be calculated by integrating

$`P^{\mathrm{max}}(L)={\displaystyle \underset{0}{\overset{\mathrm{}}{}}}xL\left(1e^{ax}\right)^{L1}e^{ax}a𝑑x`$
$`=\underset{N\mathrm{}}{lim}\{x\left(1e^{ax}\right)^L|_0^N{\displaystyle \underset{0}{\overset{N}{}}}\left(1e^{ax}\right)^L𝑑x\}.`$ (29)
An iterative splitting

$$(1e^{ax})^i=(1e^{ax})^{i1}(1e^{ax})$$
(30)
results in a summation of integrals of the type

$$\underset{0}{\overset{N}{}}(1e^{ax})^{i1}e^{ax}𝑑x=\frac{1}{ai}.$$
(31)
Thus, we obtain

$$P^{\mathrm{max}}(L)=\underset{i=1}{\overset{L}{}}\frac{1}{ai}.$$
(32)
A restriction of the statistics to values of the pressure $`P>1000`$ g cm s-2, performing an analogous calculation, results in

$$P^{\mathrm{max}}(L)=\underset{i=1}{\overset{L}{}}\frac{1}{ai}+1000\text{ g\hspace{0.17em}cm\hspace{0.17em}s}^2.$$
(33)
From the simulations we extracted the pressure distribution yielding – as expected – an exponential decrease. The prefactor is $`a=0.0009`$ for $`\mathrm{\Omega }=10`$ Hz and $`a=0.001`$ for $`\mathrm{\Omega }=2`$ Hz. The corresponding curves are plotted in Fig. 10 as dotted lines and reasonably describe to the measured values. Thus, we indirectly can conclude from the measurement that we have an exponential distribution, which indicates the occurrence of force chains.

Further insight in the mechanism of force chains is achieved through Fig. 11 which shows the average of the maximal force acting on a grain as a function of its position within the force chain. The average was performed for all force chains of length $`L=8`$. The top particle of a chain, corresponding to the maximal $`y`$-coordinate, was defined as position $`1`$. The monotonous increase of this maximal force indicates that, on the average, strain accumulates downwards along the force chain. In this respect a force chain is similar to a beam in a framework: The masses and momenta of all particles above a force chain are supported by this chain, with the consequence that lower particles of the chain are heavier loaded than upper ones. This trend is clearly visible in the data of Fig. 11. Furthermore, a force chain can shield neighbouring particles, i.e. these grains experience a much weaker force which can be observed in the simulation data as well.

Figure 12 shows the spatial pressure distribution for different rotation velocities. Here grey values code the absolute local pressure, with light pixels indicating low and dark pixels high values, in order to allow for a direct comparison of all three images.

The presented data are temporal averages. To extract the spatial distribution of this data, the two-dimensional space was coarse-grained. For each time step the positions of all particles were mapped to the grid and the resulting instantaneous local pressure was added to its related cell. Notice that the averaged data does not reflect the occurrence of force chains.

Surprisingly, for lower rotation velocities one finds that direct impacts with high relative velocities at the material surface do not result in high local pressure. Therefore, these impacts are not relevant for the fragmentation of grains in ball mills. Regions of high pressure can be found mainly near to the wall, deep inside the material. This observation is in good agreement with the experimental result of Rothkegel and Rolf Rothkegel:1992 ; Rolf . The extreme right figure in Fig. 12 reveals a different behavior. The maximum is near to the material surface where heavy impacts occur. The absolute values of the pressure are much lower than in the left figures. This observation relates to the existence of rather weak force chains.

Figs. 13 show the spatial distribution of the number of particles, which experience a pressure $`P_i>3000`$ g cm s-2. In these figures the particles which belong to a force chain are only considered in the bottom row whereas particles not belonging to a force chain are contained in the top row. Figs. 13 can be directly compared with the figures of Rothkegel Rothkegel:1992  (see Fig. 8); both are found to be in good agreement. The maxima of the distribution for $`\mathrm{\Omega }_\mathrm{I}`$ and $`\mathrm{\Omega }_{\mathrm{II}}`$ lie near to the wall of the mill, inside the material as observed in the experiment. Those maxima are an exclusive result of force chains. Grains not belonging to a force chain contribute only weakly to the comminution. For the higher velocity $`\mathrm{\Omega }_{\mathrm{III}}`$ direct impacts of particles at the surface are the dominant part. Nevertheless, absolute values are much lower than for $`\mathrm{\Omega }_\mathrm{I}`$ and $`\mathrm{\Omega }_{\mathrm{II}}`$. Again, this observation is in good agreement with the experimentRothkegel:1992 .

## 6 Conclusion

We have simulated the process of autogenous dry comminution in a ball mill by the method of molecular dynamics. To this end we developed a novel algorithm which accounts for fragmentation of particles. Throughout our simulations all particles, including the fragments, were modelled as spheres. To accomplish this idealization we introduced rather short-lived unphysical transient situations in which the fragments resulting from a cracked particle penetrate each other. The results of our simulations are available in the internet as animated sequences under http://summa.physik.hu-berlin.de/$``$kies/mill/bm.html

Employing our molecular dynamics algorithm we could satisfactorily reproduce experimental phenomena, e.g. the normalized comminution rate as a function of the normalized rotation speed and the spatial pressure distribution.

From the achieved numerical data we were able to explain an experimental result which was poorly understood so far Rothkegel:1992 : Comminution processes in ball mills occur deep inside the material rather than close to the surface where the relative collision velocity is maximal. We explained this effect by the formation of force chains. Particles which are members of one ore more force chains experience a significantly larger average force, thus substantially enhancing the fragmentation probability, as compared to particles which do not belong to any force chain.

From this we concluded that the experimentally observed and numerically achieved spatial pressure and fragmentation distributions are crucially determined by the formation of force chains. Direct collisions of particles are of minor importance for comminution. Without the existence of force chains ball mills would work much less efficient. In the context of a ball mill efficiency optimization it would be intriguing to investigate the possibility of enhancing the formation of force chains, e.g. by modifications of the geometrical properties of mills, additional application of vibration etc.

The method developed in this article may be elaborated in many different aspects:

* The numerical investigations were done in two dimensions. Ball mills are three-dimensional devices.
* Axial transport of the material was not simulated.
* All particles, including fragments, were spherical.
* Only a single-level process was simulated. Modern ball mills frequently operate in multi-level mode.
* Adhesive forces were not considered in the force law.
* Variation of material properties which originates from an increase of temperature due to comminution was not considered.

We expect that molecular dynamics will rapidly develop towards an important tool for the powder technology machinery. Due to the rapid progress in hardware and software technology, including parallel algorithms, it will be possible soon to construct and optimize powder technology facilities.

## Acknowledgement

The authors thank S. Bernotat, R. Drögemeier, E. Gommeren, G. Gudehus, H. J. Herrmann, L. Rolf, J. Schwedes and R. Weichert for discussion and for providing relevant literature. The work was supported by German Science Foundation through project Po472/3-2.


text:
# On a Factorization of Symmetric Matrices and Antilinear Symmetries

## 1 Introduction

One of the interesting results of linear algebra is that every square matrix may be factored into the product of two symmetric matrices . The factorization of symmetric matrices into a product of a square matrix and its transpose is however less known. In fact, there seems to be no mention of this factorization in modern texts on linear algebra. The purpose of this note is to present a simple derivation of this particular factorization of symmetric matrices and to discuss its application in establishing a uniqueness property of certain antilinear operators.

## 2 Notation and Definitions

In this note we shall express the complex-conjugate, the transpose, and the conjugate-transpose (adjoint) of a matrix (an operator) $`m`$ by $`\overline{m}`$, $`m^T`$, and $`m^{}`$, respectively, and identify the elements of $`^{\mathrm{}}`$ by columns of $`\mathrm{}`$ complex numbers. Then the Euclidean inner-product of two vectors $`\stackrel{}{w}_1`$ and $`\stackrel{}{w}_2`$ takes the form $`\stackrel{}{w}_2^{}\stackrel{}{w}_1`$.

* Definition 1: Let $`S`$ be a set, then the Kronecker delta function $`\delta :S^2\{0,1\}`$ is defined by

$$a,bS,\delta (a,b)=\delta _{ab}:=\{\begin{array}{ccc}1& \mathrm{for}& a=b\\ 0& \mathrm{for}& ab.\end{array}$$
* Definition 2: A function $`𝒳:`$ acting in a complex vector space $``$ is said to be an antilinear operator if for all $`x,y`$ and $`\varphi ,\psi `$, $`𝒳(x\varphi +y\psi )=\overline{x}𝒳\varphi +\overline{y}𝒳\psi `$, where $`𝒳\varphi `$ means $`𝒳(\varphi )`$.
* Definition 3: An antilinear operator $`𝒳:`$ acting in a complex inner-product space $``$ with inner-product $`(,)`$ is said to be symmetric or Hermitian if for all $`\varphi ,\psi `$, $`(𝒳\varphi ,\psi )=(𝒳\psi ,\varphi )`$.
* Definition 4: A linear operator $`H:`$ acting in a complex inner-product space $``$ is said to have a symmetry generated by a function $`X:`$ or simply a $`X`$-symmetry if $`H`$ and $`X`$ commute, i.e., $`[H,X]:=HXXH=0`$, where $`0`$ stands for the zero operator. A symmetry generated by an antilinear operator is called an antilinear symmetry.
* Definition 5: Let $`H:`$ be a linear operator acting in a complex inner-product space $``$ and $`G:`$ be an Hermitian invertible linear or antilinear operator. Then $`H`$ is said to be $`G`$-Hermitian  or $`G`$-pseudo-Hermitian  if $`H^{}=GHG^1`$.

## 3 Motivation: Consequences of antilinear symmetries

Consider a diagonalizable linear operator $`H:`$ acting in a finite-dimensional complex inner-product space $``$ with inner-product $`(,)`$. Let $`n`$ label the eigenvalues $`E_n`$ of $`H`$, $`\mu _n`$ be the multiplicity of $`E_n`$, and $`\psi _{n,a}`$ be the eigenvectors corresponding to the eigenvalue $`E_n`$ where $`a\{1,2,\mathrm{},\mu _n\}`$ is the degeneracy label. Then it is well-known  that the adjoint $`H^{}`$ of $`H`$ is diagonalizable; the eigenvalues $`\stackrel{~}{E}_n`$ of $`H^{}`$ are complex conjugate of those of $`H`$, i.e., $`\stackrel{~}{E}_n=\overline{E}_n`$; the multiplicity of $`\stackrel{~}{E}_n`$ is equal to $`\mu _n`$; and one can choose the eigenvectors $`\varphi _{n,a}`$ of $`H^{}`$ in such a way that for all spectral labels $`m,n`$ and degeneracy labels $`a,b`$, $`(\varphi _{n,a},\psi _{n,a})=\delta _{n,m}\delta _{a,b}`$. Clearly, both sets of eigenvectors $`\psi _{n,a}`$ of $`H`$ and eigenvectors $`\varphi _{a,n}`$ of $`H^{}`$ form bases of $``$; $`\{\psi _{n,a},\varphi _{n,a}\}`$ is a complete biorthonormal system. Recently , we have shown that if the eigenvalues of $`H`$ are real, then $`H`$ has an antilinear symmetry. More generally, we proved the following theorem.

* Theorem 1: The presence of an antilinear symmetry of $`H`$ is a necessary and sufficient condition for the eigenvalues of $`H`$ to either be real or come in complex conjugate pairs.

The proof of Theorem 1 uses the following lemma.

* Lemma 1: Every diagonalizable linear operator $`H:`$ acting in a finite-dimensional complex inner-product space $``$ is $`𝒯`$-Hermitian,

$$H^{}=𝒯H𝒯^1,$$
(1)
for some Hermitian, invertible, antilinear operator $`𝒯:`$.

It turns out  that any such $`𝒯`$ may be expressed in terms of the eigenvectors $`\varphi _{n,a}`$ of $`H^{}`$ according to

$$\zeta ,𝒯\zeta =\underset{n}{}\underset{a,b=1}{\overset{\mu _n}{}}c_{ba}^{(n)}(\varphi _{n,a},\zeta )\varphi _{n,b},$$
(2)
where $`c_{ab}^{(n)}`$ are the entries of symmetric invertible $`\mu _n\times \mu _n`$ matrices $`c^{(n)}`$.

Note that Theorem 1 and Lemma 1 have infinite-dimensional generalizations for linear operators $`H`$ admitting a complete biorthonormal system of eigenvectors .

Next, consider a general basis transformation,

$$\varphi _{n,a}\varphi _{n,a}^{}:=\underset{b=1}{\overset{\mu _n}{}}v_{ba}^{(n)}\varphi _{n,b}$$
(3)
where $`v_{ab}^{(n)}`$ are the entries of an invertible $`\mu _n\times \mu _n`$ matrix $`v^{(n)}`$. In terms of the transformed basis vectors $`𝒯`$ has the form: $`\zeta `$, $`𝒯\zeta =_n_{a,b=1}^{\mu _n}c_{ba}^{{}_{}{}^{}(n)}(\varphi _{n,a}^{},\zeta )\varphi _{n,b}^{}`$, where $`c^{{}_{}{}^{}(n)}`$ are related to $`c^{(n)}`$ according to $`c^{(n)}=v^{(n)}c^{{}_{}{}^{}(n)}v^{(n)T}.`$ This equation indicates that the issue of the uniqueness of $`𝒯`$ for a given $`H`$ is related to whether one can find for each $`c^{(n)}`$ an invertible matrix $`v^{(n)}`$ such that $`c^{(n)}=v^{(n)}v^{(n)T}`$. In the remainder of this note we shall give a proof of the fact that this is indeed possible, and one can transform to a basis in which $`𝒯`$ has the (canonical) form (2) with $`c_{ab}^{(n)}=\delta _{ab}`$ for all $`n`$.

## 4 Factorization of Symmetric Matrices

* Theorem 2: A square matrix $`c`$ is symmetric if and only if it can be written as $`c=vv^T`$ for some square matrix $`v`$.
* Proof: If $`c=vv^T`$ then clearly $`c`$ is symmetric. To prove the converse we use induction on the dimension $`n`$ of the matrix $`c`$. For $`n=1`$, $`c=vv^T=v^2`$ is trivially satisfied by letting $`v:=\sqrt{c}`$. By induction hypothesis we assume that for all $`k\{2,\mathrm{},n\}`$, every $`k\times k`$ symmetric matrix $`c`$ can be written in the form $`c=vv^T`$ for some $`k\times k`$ matrix $`v`$. Now let $`C`$ be an $`(n+1)\times (n+1)`$ symmetric matrix. Then $`C`$ has at least one eigenvector  i.e., there are $`\stackrel{}{e}^{n+1}\{\stackrel{}{0}\}`$ and $`\lambda `$ such that

$$C\stackrel{}{e}=\lambda \stackrel{}{e}.$$
(4)
Now let $`𝒱:=\{\stackrel{}{w}^{n+1}|\stackrel{}{w}^{}\overline{\stackrel{}{e}}=0\}`$ be the orthogonal complement of $`\overline{\stackrel{}{e}}`$. Clearly $`𝒱`$ is an $`n`$-dimensional vector subspace of $`^{n+1}`$. Next, consider the following two possibilities.

+ $`\stackrel{}{e}𝒱`$. In this case, choose a basis $`\{\stackrel{}{e}_1,\stackrel{}{e}_n,\mathrm{},\stackrel{}{e}_n\}`$ of $`𝒱`$ and let $`\stackrel{}{e}_{n+1}:=\stackrel{}{e}`$. Then $`\{\stackrel{}{e}_1,\stackrel{}{e}_2,\mathrm{},\stackrel{}{e}_n,\stackrel{}{e}_{n+1}\}`$ is a basis of $`^{n+1}`$ and the matrix $`A:=(\stackrel{}{e}_1,\stackrel{}{e}_n,\mathrm{},\stackrel{}{e}_n,\stackrel{}{e}_{n+1})`$ is invertible. Note that for all $`\mathrm{}\{1,2,\mathrm{},n\}`$, $`\stackrel{}{e}_{\mathrm{}}V`$, and $`\stackrel{}{e}_{\mathrm{}}^T\stackrel{}{e}_{n+1}=0`$. This in turn implies that the matrix $`A^TCA`$ which is symmetric has the block form

    $$A^TCA=\left(\begin{array}{cc}\stackrel{~}{c}& \stackrel{}{0}\\ \stackrel{}{0}^T& \lambda ^2\end{array}\right),$$
    (5)
    where $`\stackrel{~}{c}`$ is a symmetric $`n\times n`$ matrix. By induction hypothesis there is an $`n\times n`$ matrix $`\stackrel{~}{v}`$ such that $`\stackrel{~}{c}=\stackrel{~}{v}\stackrel{~}{v}^T`$. Now let $`B`$ be the $`(n+1)\times (n+1)`$ matrix

    $$B:=\left(\begin{array}{cc}\stackrel{~}{v}^T& \stackrel{}{0}\\ \stackrel{}{0}^T& \lambda \end{array}\right),$$
    (6)
    and $`V:=(BA^1)^T`$. Then in view of (5) and (6), $`B^TB=A^TCA`$ and

    $$VV^T=(BA^1)^TBA^1=A^{1T}B^TBA^1=C.$$

    This completes the proof for case (i).
    + $`\stackrel{}{e}𝒱`$, i.e., $`\stackrel{}{e}^T\stackrel{}{e}=0`$. In this case, let $`𝒱^{}:=\{\stackrel{}{w}𝒱|\stackrel{}{w}^{}\stackrel{}{e}=0\}`$ be the orthogonal complement of $`\stackrel{}{e}`$ in $`𝒱`$, $`\{\stackrel{}{e}_1,\stackrel{}{e}_2,\mathrm{},\stackrel{}{e}_{n1}\}`$ be a basis of $`𝒱^{}`$, $`\stackrel{}{e}_n=\overline{\stackrel{}{e}}`$, and $`\stackrel{}{e}_{n+1}=\stackrel{}{e}`$. Then $`\{\stackrel{}{e}_1,\stackrel{}{e}_2,\mathrm{},\stackrel{}{e}_n,\stackrel{}{e}_{n+1}\}`$ is a basis of $`^{n+1}`$ and the matrix $`A^{}:=(\stackrel{}{e}_1,\stackrel{}{e}_n,\mathrm{},\stackrel{}{e}_n,\stackrel{}{e}_{n+1})`$ is invertible. Note that for all $`\mathrm{}\{1,2,\mathrm{},n1\}`$, $`\stackrel{}{e}_{\mathrm{}}V`$, and $`\stackrel{}{e}_{\mathrm{}}^T\stackrel{}{e}_{n+1}=0`$. Furthermore, $`\stackrel{}{e}_{n+1}^T\stackrel{}{e}_{n+1}=\stackrel{}{e}^T\stackrel{}{e}=0`$ and $`\alpha :=\stackrel{}{e}_n^T\stackrel{}{e}_{n+1}=\stackrel{}{e}^{}\stackrel{}{e}^+`$. In view of these relations and (4),

    $$C^{}:=A^{}_{}{}^{}TCA^{}=\left(\begin{array}{cccccc}\stackrel{~}{c}_{1,1}^{}& \stackrel{~}{c}_{1,2}^{}& \mathrm{}& \stackrel{~}{c}_{1,n1}^{}& \stackrel{~}{c}_{1,n}^{}& 0\\ \stackrel{~}{c}_{2,1}^{}& \stackrel{~}{c}_{2,2}^{}& \mathrm{}& \stackrel{~}{c}_{2,n1}^{}& \stackrel{~}{c}_{2,n}^{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& 0\\ \stackrel{~}{c}_{n1,1}^{}& \stackrel{~}{c}_{n1,2}^{}& \mathrm{}& \stackrel{~}{c}_{n1,n1}^{}& \stackrel{~}{c}_{n1,n}^{}& 0\\ \stackrel{~}{c}_{n,1}^{}& \stackrel{~}{c}_{n,2}^{}& \mathrm{}& \stackrel{~}{c}_{n,n1}^{}& \stackrel{~}{c}_{n,n}^{}& \lambda \alpha \\ 0& 0& \mathrm{}& 0& \lambda \alpha & 0\end{array}\right),$$
    (7)
    where $`\stackrel{~}{c}_{i,j}^{}:=\stackrel{}{e}_i^TC\stackrel{}{e}_j`$ are the entries of a symmetric $`n\times n`$ matrix $`\stackrel{~}{c}^{}`$. Now if $`\lambda =0`$, $`C^{}`$ is block-diagonal and the argument given in case (i) leads to a proof of the theorem. This leaves the case $`\lambda 0`$. In this case, let $`A:=A^{}D`$ where $`D`$ is the $`(n+1)\times (n+1)`$ matrix

    $$D=\left(\begin{array}{ccccccc}1& 0& 0& \mathrm{}& 0& 0& x_1\\ 0& 1& 0& \mathrm{}& 0& 0& x_2\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& 0& \mathrm{}& 1& 0& x_{n1}\\ 0& 0& 0& \mathrm{}& 0& 1& x_n\\ y_1& y_2& y_3& \mathrm{}& y_{n1}& y_n& 0\end{array}\right),$$
    (8)
    with $`x_1,x_2,\mathrm{}x_{n1}`$ being arbitrary complex numbers,

    $$x_n:=(\lambda \alpha )^1,$$
    (9)
    and for all $`i\{1,2,\mathrm{},n\}`$

    $$y_i:=\underset{j=1}{\overset{n}{}}\stackrel{~}{c}_{i,j}^{}x_j.$$
    (10)
    Then a simple computation shows that

    $$A^TCA=D^TC^{}D=\left(\begin{array}{cc}\stackrel{~}{c}_{ij}& \stackrel{}{0}\\ \stackrel{}{0}^T& \lambda ^{}_{}{}^{}2\end{array}\right),$$

    where $`\stackrel{~}{c}_{i,j}`$ are the entries of a symmetric $`n\times n`$ matrix $`\stackrel{~}{c}`$ and $`\lambda ^{}`$. Therefore, we can use the argument given in case (i) to show the existence of an $`(n+1)\times (n+1)`$ matrix $`B`$ satisfying

    $$A^TCA=B^TB.$$
    (11)
    The proof of the theorem will be complete if we show that the matrix $`A=A^{}D`$ is invertible. Because $`A^{}`$ is invertible, it suffices to show the existence of $`x_1,x_2,\mathrm{}x_{n1}`$ for which $`detD0`$. We can use the properties of the determinant and Equations (9) and (10) to compute

    $`detD`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}x_iy_i`$ (12)
    $`=`$ $`{\displaystyle \underset{i=1}{\overset{n1}{}}}\stackrel{~}{c}_{i,i}^{}x_i^22{\displaystyle \underset{i<j=1}{\overset{n1}{}}}\stackrel{~}{c}_{i,j}^{}x_ix_j+2(\lambda \alpha )^1{\displaystyle \underset{i=1}{\overset{n1}{}}}\stackrel{~}{c}_{n,i}^{}x_i(\lambda \alpha )^2\stackrel{~}{c}_{n,n}^{}.`$
    Suppose that for all values of $`x_1,x_2,\mathrm{}x_{n1}`$, $`detD=0`$. This implies that $`\stackrel{~}{c}^{}=0`$ in which case $`C^{}`$ will have the form

    $$C^{}=\left(\begin{array}{cccccc}0& 0& \mathrm{}& 0& 0& 0\\ 0& 0& \mathrm{}& 0& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 0& 0& 0\\ 0& 0& \mathrm{}& 0& 0& \lambda \alpha \\ 0& 0& \mathrm{}& 0& \lambda \alpha & 0\end{array}\right).$$

    By induction hypothesis we have a $`2\times 2`$ matrix $`m`$ satisfying

    $$\left(\begin{array}{cc}0& \lambda \alpha \\ \lambda \alpha & 0\end{array}\right)=m^Tm.$$

    Therefore, setting

    $$B:=\left(\begin{array}{cccccc}0& 0& \mathrm{}& 0& 0& 0\\ 0& 0& \mathrm{}& 0& 0& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 0& 0& 0\\ 0& 0& \mathrm{}& 0& m_{1,2}& m_{1,2}\\ 0& 0& \mathrm{}& 0& m_{2,1}& m_{2,2}\end{array}\right)$$

    and $`V:=(BA^{}_{}{}^{}1)^T`$, we have $`C^{}=B^TB`$ and $`C=VV^T`$. If $`\stackrel{~}{c}^{}0`$ then there are values of $`x_1,x_2,\mathrm{}x_{n1}`$ for which $`D`$ is a nonsingular matrix and $`A=A^{}D`$ is invertible. Hence we can set $`V:=(BA^1)^T`$ and use (11) to show that $`C=VV^T`$. $`\mathrm{}`$

## 5 Concluding remarks

1. The factorization $`c=vv^T`$ established in Theorem 2 is invariant under the transformation $`vv^{}=vo`$ where $`o`$ is an arbitrary (complex) orthogonal matrix. In particular, one may choose $`o`$ so that the factorizing matrix $`v^{}`$ has a simple form.
2. In view of the discussion of Section 3, one has the following consequence of Theorem 2.

* Corollary 1: Up to basis transformations (3), there is a unique antilinear operator $`𝒯`$ satisfying Equation (1), namely

    $$\zeta ,𝒯\zeta =\underset{n}{}\underset{a=1}{\overset{\mu _n}{}}(\varphi _{n,a},\zeta )\varphi _{n,a},$$
    (13)
3. For a self-adjoint linear operator $`H`$, one can set $`\varphi _{n,a}=\psi _{n,a}`$ and use the completeness of the eigenvectors $`\psi _{n,a}`$ and (13) to deduce $`𝒯^2=I`$, where $`I`$ is the identity operator. Furthermore, noting that in this case (1) is equivalent to $`𝒯`$-symmetry of $`H`$, one can prove the following.

* Corollary 2: Every self-adjoint linear operator $`H`$ has an antilinear symmetry generated by a Hermitian, invertible, antilinear operator $`𝒯`$ satisfying $`𝒯^2=I`$.

## Acknowledgment

This project was supported by the Young Researcher Award Program (GEBIP) of the Turkish Academy of Sciences.