text:
# Phonons in a nanoparticle mechanically coupled to a substrate

## I Introduction

There is currently great interest in properties of nanometer-scale mechanical systems, such as cantilevers, nanoparticles, and resonators.roukes  Because of the extremely small size and volume-to-surface ratio of these systems, interactions with their surroundings can dramatically alter their properties. In particular, it is well known that the vibrational spectrum of an isolated nanometer-scale crystal, being discrete, is qualitatively different than that of the same bulk material, leading to important changes in any property dependent on the phonon density of states (DOS). The differences between the vibrational DOS in a nanoparticle and a bulk solid are most evident at low frequencies: A spherical nanoparticle with diameter $`d`$ and characteristic bulk sound velocity $`v`$ can not support a mode with frequency less than about $`2\pi v/d`$. Thus, an acoustic “gap” in the low-energy phonon spectrum is present in contrast with that of the bulk, which has a continuous spectrum down to zero energy. However, mechanical interaction with the environment will modify the discrete nature of the modes.

In an interesting experiment by Yang et al.,yang  the phonon DOS deep inside this gap was measured in insulating $`\mathrm{Y}_2\mathrm{O}_3`$ nanoparticles. The experiment consisted of nanoparticles whose sizes ranged from 7 to 23 nm in diameter and was performed by measuring the nonradiative lifetimes of an excited electronic state of a $`\mathrm{Eu}^{3+}`$ dopant. The lowest supported mode or Lamb mode for nanoparticles of these sizes is about $`10\mathrm{cm}^1`$. At $`3\mathrm{cm}^1`$ the DOS measured was more than 100 times smaller than that of bulk $`\mathrm{Y}_2\mathrm{O}_3`$ (at 3 cm<sup>-1</sup>).

In this paper we propose and investigate a mechanism that could be responsible for the observed broadening of the nanoparticles’ phonon modes. Several possible broadening mechanisms could be responsible for the observed effect. For example, anharmonicity leads to broadening and, therefore, to a low-energy DOS, but anharmonicity is ineffective at low energy and was found to be too small to account for the experiment.vadim  Another possibility could be adsorbed “dirt” on the outside of the nanoparticle. This might lower the Q factor of the nanoparticle, regarding it as a resonator, reflecting a broadening of the vibrational modes. A third, and in our opinion more likely mechanism, follows from the realization that these nanoparticles are not isolated, but instead were prepared in a powered form. Thus, each nanoparticle is in weak contact with a cluster of other nanoparticles. Because the cluster is relatively large, including several hundred nanoparticles, each nanoparticle is mechanically coupled to a reservoir that has a continuous vibrational DOS at low energy. This interaction broadens the modes and allows phonons in the nanoparticle to escape and be absorbed into the cluster. We investigate the affect this mechanical environmental interaction has on the nanoparticle’s phonon DOS.

Because we are only interested in determining the correct origin of broadening, and do not hope to be able to exactly reproduce the experimental results of Ref. yang, , we propose the following simplified model: The cluster of nanoparticles is replaced by a semi-infinite elastic substrate, and one nanoparticle is placed in weak mechanical contact with it. The weak contact is imagined to be a few atomic bonds or small neck of material, which we model by a harmonic spring. For simplicity, we take the substrate and the nanoparticle to be made out of the same isotropic elastic material. Because we are interested in the low energy regime, continuum elasticity theory will be used to describe the nanoparticle and the substrate. After defining and analyzing our simplified model, in Sections III through V, we explain in Section VI how the model can be adapted to address the experiment of Ref. yang, , and good agreement is obtained.

The simple model we study is related to, but different than, models used to study energy relaxation by molecules adsorbed on surfaces.zangwill  However, in surface science the interest is usually in the relaxation of rigid translational motion, rotational motion, or simple internal vibrations of adsorbates. In contrast, we investigate the broadening of complex internal vibrational modes of much larger objects (which are crystalline). Also, our work has much in common with that of Gurevich and Schober,g and s  where many of the same considerations and modeling were used to study the Lamb-mode decay rate, of nanoparticles caused by both anharmonicity and coupling to an enviroment of other nanoparticles.

## II NANOPARTICLE AND SUBSTRATE MODEL

As mentioned in the introduction, the model we study is that of a single nanoparticle in weak mechanical contact with a semi-infinite substrate. Linear elasticity theory will be used to describe the phonons of this system. We assume the nanoparticle and substrate to be made of an isotropic non-polar material. Because we take the nanoparticle and substrate to be made of the same material, we will use the same density $`\rho `$ and Lamé coefficients $`\lambda `$ and $`\mu `$ for both. The Lagrangian for the entire system is given by

$$L=_𝖵d^3r\left[\frac{1}{2}\rho (_t𝐮)^2\frac{1}{2}\lambda u_{ii}^2\mu u_{ij}^2\right],$$
(1)
where $`𝐮(𝐫,t)`$ is the displacement field, and

$$u_{ij}(_iu_j+_ju_i)/2$$
(2)
is the strain tensor. $`𝖵`$ is the volume of the nanoparticle, substrate, and connecting material, as shown in Fig. 1.

Because the Lagrangian density is local, the integration volume in (1) can be split into three independent parts: the nanoparticle, the substrate, and the connecting region. In the limit of weak coupling (diameter $`d_\mathrm{c}`$ of connecting region much smaller than $`d`$), the surface area on the nanoparticle and substrate over which the actual boundary conditions differs from stress-free conditions are negligible, and the Hamiltonian can be written as

$$H=H_{\mathrm{nano}}+H_{\mathrm{sub}}+\delta H,$$
(3)
where $`H_{\mathrm{nano}}`$ is the Hamiltonian for an isolated nanoparticle (with stress-free boundaries), $`H_{\mathrm{sub}}`$ is that for an isolated substrate, and $`\delta H`$ is the interaction between the two. The connecting region is taken to be a few atomic bonds or small neck of material, as shown in Fig. 1. We further approximate this mechanical coupling to be a harmonic spring potential

$$\delta H=\frac{1}{2}K:\left[u_{\mathrm{nano}}^z(𝐫_0)u_{\mathrm{sub}}^z(𝐫_0)\right]^2:,$$
(4)
where $`K`$ is an effective spring constant, and $`u_{\mathrm{nano}}^z`$ and $`u_{\mathrm{sub}}^z`$ are the $`z`$ components of the displacement field of the nanoparticle and substrate at the point of contact, $`𝐫_0`$. We take the $`z`$ direction to be along the upward pointing normal to the substrate surface. The Hamiltonians we will introduce below for $`H_{\mathrm{nano}}`$ and $`H_{\mathrm{sub}}`$ are normal-ordered; therefore, it is necessary to normal order $`\delta H`$ as well. This operation is denoted by the colons in Eq. (4).

Our analysis will require the vibrational normal modes and spectra of the isolated nanoparticle and semi-infinite substrate, calculated with stress-free boundary conditions. The long-wavelength modes of interest here may be obtained from elasticity theory, to which we now turn.

### II.1 Isolated Nanoparticle

Here we derive the normal modes of an isolated elastic sphere. The method we shall use is different than (but equivalent to) that used in the classic paper by Lamb,lamb  but is better suited for our purposes. The equation of motion of (1) is

$$_t^2𝐮v_\mathrm{l}^2\mathbf{}(\mathbf{}𝐮)+v_\mathrm{t}^2\mathbf{}\times \mathbf{}\times 𝐮=0,$$
(5)
where $`v_\mathrm{l}\sqrt{(\lambda +2\mu )/\rho }`$ is the bulk longitudinal sound velocity and $`v_\mathrm{t}\sqrt{\mu /\rho }`$ is the transverse velocity. To solve Eq. (5) the displacement field can be decomposed into longitudinal and transverse parts,

$$𝐮=𝐮_\mathrm{l}+𝐮_\mathrm{t},$$
(6)
where

$$\mathbf{}\times 𝐮_\mathrm{l}=0$$
(7)
and

$$\mathbf{}𝐮_\mathrm{t}=0.$$
(8)
With harmonic time dependence, the equation of motion Eq. (5) then separates into two vector Helmholtz equations for the longitudinal and transverse parts,

$$\left(^2+p^2\right)𝐮_\mathrm{l}=0,p\omega /v_\mathrm{l}$$
(9)
and

$$\left(^2+q^2\right)𝐮_\mathrm{t}=0,q\omega /v_\mathrm{t}.$$
(10)
The longitudinal equation (9) can be solved by introducing a scalar potential

$$𝐮_\mathrm{l}=\mathbf{}\varphi ^{(p)},$$
(11)
where $`\varphi ^{(p)}`$ is a solution of the scalar Helmholtz equation $`\left(^2+p^2\right)\varphi ^{(p)}=0`$. The transverse equation (10) has two linearly independent solutions, $`𝐮_\mathrm{t}=𝐌\mathrm{and}𝐍,`$ where

$$𝐌=\mathbf{}\varphi ^{(q)}\times 𝐫$$
(12)
and

$$𝐍=\frac{1}{q}\mathbf{}\times 𝐌.$$
(13)
Here $`\varphi ^{(q)}`$ is a solution of $`\left(^2+q^2\right)\varphi ^{(q)}=0`$. The prefactor $`1/q`$ is included for dimensional convenience. The scalar Helmholtz equations are separable in spherical coordinates and the solutions can be written as

$$\varphi _{lm}(𝐫)j_l(kr)Y_{lm}(\theta ,\phi ),k=p,q$$
(14)
where

$$j_l(x)\sqrt{\frac{\pi }{2x}}J_{l+\frac{1}{2}}(x)$$
(15)
is a spherical Bessel function of the first kind (regular at origin) and

$$Y_{lm}(\theta ,\phi )(1)^m\sqrt{\frac{2l+1}{4\pi }\frac{(lm)!}{(l+m)!}}P_{lm}(\mathrm{cos}\theta )e^{im\phi }.$$
(16)
Here

$$P_{lm}(x)(1x^2)^{\frac{m}{2}}\frac{^m}{x^m}P_l(x),$$
(17)
where $`P_l(x)`$ are Legendre polynomials.

Now we use the $`\varphi _{lm}`$ to construct three linearly independent solutions of (5),

$`𝐋_{lm}`$ $``$ $`\frac{1}{p}\mathbf{}\varphi _{lm}(pr),`$ (18)
$`𝐌_{lm}`$ $``$ $`\mathbf{}\varphi _{lm}(qr)\times 𝐫,`$ (19)
$`𝐍_{lm}`$ $``$ $`\frac{1}{q}\mathbf{}\times 𝐌_{lm}.`$ (20)
The general solution is a linear combination of $`𝐋_{lm}`$, $`𝐌_{lm}`$, and $`𝐍_{lm}`$,

$$𝐮(𝐫)=\underset{lm}{}\left[a_{lm}𝐋_{lm}+b_{lm}𝐌_{lm}+c_{lm}𝐍_{lm}\right].$$
(21)
Although they are linearly independent, the vector fields $`𝐋_{lm}`$, $`𝐌_{lm}`$, and $`𝐍_{lm}`$ are not orthogonal in space. However, they can be written in terms of orthogonal vector spherical harmonics $`𝐏_{lm},𝐁_{lm}`$, and $`𝐂_{lm}`$, defined as

$$𝐏_{lm}(\mathrm{\Omega })Y_{lm}(\mathrm{\Omega })𝐞_r$$
(22)
$$𝐁_{lm}(\mathrm{\Omega })\frac{1}{\sqrt{l(l+1)}}\left(_\theta Y_{lm}(\mathrm{\Omega })𝐞_\theta +\frac{imY_{lm}(\mathrm{\Omega })}{\mathrm{sin}\theta }𝐞_\phi \right)$$
(23)
$$𝐂_{lm}(\mathrm{\Omega })\frac{1}{\sqrt{l(l+1)}}\left(\frac{imY_{lm}(\mathrm{\Omega })}{\mathrm{sin}\theta }𝐞_\theta _\theta Y_{lm}(\mathrm{\Omega })𝐞_\phi \right),$$
(24)
with the following properties

$$𝑑\mathrm{\Omega }𝐗_{lm}^{}𝐗_{l^{}m^{}}=\delta _{ll^{}}\delta _{mm^{}}$$
(25)
for $`𝐗𝐁,𝐂,𝐏`$ and

$$𝑑\mathrm{\Omega }𝐗_{lm}^{}𝐗_{l^{}m^{}}^{}=0$$
(26)
for $`𝐗𝐗^{}`$. Expressed in terms of vector spherical harmonics, $`𝐋_{lm},𝐌_{lm}`$ and $`𝐍_{lm}`$ are given by

$$𝐋_{lm}=j_l^{}(pr)𝐏_{lm}(\mathrm{\Omega })+\frac{\sqrt{l(l+1)}}{pr}j_l(pr)𝐁_{lm}(\mathrm{\Omega }),$$
(27)
$`𝐌_{lm}=\sqrt{l(l+1)}j_l(qr)𝐂_{lm}(\mathrm{\Omega }),`$ (28)
and

$`𝐍_{lm}`$ $`=`$ $`{\displaystyle \frac{l(l+1)}{qr}}j_l(qr)𝐏_{lm}(\mathrm{\Omega })+{\displaystyle \frac{\sqrt{l(l+1)}}{qr}}`$ (29)
$`\times `$ $`\left[j_l(qr)+qrj_l^{}(qr)\right]𝐁_{lm}(\mathrm{\Omega }),`$
where prime denotes differentiation with respect to the argument.

Next we impose stress-free boundary conditions

$$\sigma _{ij}n_j=0$$
(30)
at the surface $`r=R`$ of the nanoparticle. Here $`𝐧`$ is an outward pointing normal vector and $`\sigma _{ij}`$ is the strain tensor. In an isotropic elastic continuum,

$$\sigma _{ij}=\lambda (\mathbf{}𝐮)\delta _{ij}+2\mu u_{ij}.$$
(31)
In spherical coordinates (30) implies

$$\sigma _{rr}=\sigma _{\theta r}=\sigma _{\phi r}=0.$$
(32)
The three conditions (32) require that

$$\lambda (\mathbf{}𝐮)+2\mu u_{rr}=0,$$
(33)
$$u_{\theta r}=0,$$
(34)
and

$$u_{\phi r}=0.$$
(35)
In terms of the displacement field,

$`u_{rr}`$ $`=`$ $`_ru_r,`$ (36)
$`u_{\theta r}`$ $`=`$ $`\frac{1}{2}(_ru_\theta \frac{1}{r}u_\theta +\frac{1}{r}_\theta u_r),`$ (37)
$`u_{\phi r}`$ $`=`$ $`\frac{1}{2}(\frac{1}{r\mathrm{sin}\theta }_\phi u_r+_ru_\phi \frac{1}{r}u_\phi ).`$ (38)
The boundary condition equations (32) then become

$`a_{lm}\left[\lambda pj_l(pR)Y_{lm}+2\mu pj_l^{\prime \prime }(pR)Y_{lm}\right]`$
$`+c_{lm}2\mu l(l+1)𝖤Y_{lm}=0,`$ (39)
$`a_{lm}2𝖣_\theta Y_{lm}+b_{lm}im𝖤Y_{lm}\mathrm{csc}\theta `$
$`+c_{lm}𝖥_\theta Y_{lm}=0,`$ (40)
and

$`a_{lm}2im𝖣Y_{lm}\mathrm{csc}\theta +b_{lm}𝖤_\theta Y_{lm}`$
$`+c_{lm}im𝖥Y_{lm}\mathrm{csc}\theta =0,`$ (41)
where

$`𝖣`$ $`{\displaystyle \frac{j_l^{}(pR)}{R}}{\displaystyle \frac{j_l(pR)}{pR^2}},`$ (42)
$`𝖤`$ $`{\displaystyle \frac{j_l^{}(qR)}{R}}{\displaystyle \frac{j_l(qR)}{qR^2}},`$ (43)
$`𝖥`$ $`qj_l^{\prime \prime }(qR)+{\displaystyle \frac{l(l+1)}{qR^2}}j_l(qR){\displaystyle \frac{2j_l(qR)}{qR^2}}.`$ (44)
Finally, we rewrite (39) though (41) in matrix form as

$$\left(\begin{array}{ccc}\lambda pj_l(pR)Y_{lm}+2\mu pj_l^{\prime \prime }(pR)Y_{lm}& 0& 2\mu l(l+1)𝖤Y_{lm}\\ 2𝖣_\theta Y_{lm}& im𝖤Y_{lm}\mathrm{csc}\theta & 𝖥_\theta Y_{lm}\\ 2im𝖣Y_{lm}\mathrm{csc}\theta & 𝖤_\theta Y_{lm}& im𝖥Y_{lm}\mathrm{csc}\theta \end{array}\right)\left(\begin{array}{c}a_{lm}\\ b_{lm}\\ c_{lm}\end{array}\right)=0.$$

For a nontrivial solution of (39) through (41) to exist, the determinant of the above matrix must vanish. Taking the determinant and simplifying we find

$$\left[\lambda pj_l(pR)+2\mu pj_l^{\prime \prime }(pR)\right]𝖤𝖥4\mu l(l+1)𝖤^2𝖣=0.$$
(45)
This implies that either

$$𝖤=\frac{j_l^{}(qR)}{R}\frac{j_l(qR)}{qR^2}=0$$
(46)
or

$$\left[\lambda pj_l(pR)+2\mu pj_l^{\prime \prime }(pR)\right]𝖥4\mu l(l+1)𝖤𝖣=0.$$
(47)
If (46) is met, then this imposes certain constraints on $`a_{lm},b_{lm}`$ and $`c_{lm}`$, which require $`a_{lm}=c_{lm}=0`$. This can easily be seen in the above matrix by setting $`𝖤=0`$. If (47) is met, $`b_{lm}`$ has to be zero. In conclusion, we have two branches of vibrational modes. The branch in which (46) is satisfied,

$$𝐮(𝐫)=b_{lmn}𝐌_{lmn}(𝐫),$$
(48)
are referred to as the torsional modes, where $`n`$ specifies the radial quantum number. The other branch is found when (47) is satisfied,

$$𝐮(𝐫)=a_{lmn}𝐋_{lmn}(𝐫)+c_{lmn}𝐍_{lmn}(𝐫),$$
(49)
which are the spheroidal modes.

To quantize the vibrational modes we write the displacement field asfootnote1

$$𝐮_{\mathrm{nano}}(𝐫)=\underset{J}{}\sqrt{\frac{\mathrm{}}{2\rho \omega _J}}\left[a_J𝚿_J(𝐫)+a_J^{}𝚿_J^{}(𝐫)\right],$$
(50)
where

$$J=[\mathrm{S}\mathrm{or}\mathrm{T},n,l,m]$$
(51)
is a label uniquely specifying a nanoparticle eigenmode. The first entry S or T specifies whether the mode is in the spheroidal or torsional branch, respectively. $`n`$ is the radial quantum number and $`l`$ and $`m`$ are the usual angular momentum quantum numbers. $`a`$ and $`a^{}`$ are phonon annihilation and creation operators which satisfy the Bose commutation relation

$$[a_J,a_J^{}^{}]=\delta _{JJ^{}}.$$
(52)
The $`𝚿_J`$ are vibrational eigenvectors normalized such that

$$_V|𝚿_J(𝐫)|^2d^3r=1,$$
(53)
where $`V`$ is the volume of the nanoparticle. Assuming (without proof) that the modes $`𝚿_J`$ form a complete set,

$$\underset{J}{}\mathrm{\Psi }_J^i(𝐫)\mathrm{\Psi }_J^j(𝐫^{})=\delta ^{ij}\delta (𝐫𝐫^{}),$$
(54)
it can easily be shown that $`𝐮`$ satisfies the correct equal-time canonical commutation relation with $`𝝅\rho _t𝐮`$, namely

$$[u^i(𝐫),\pi ^j(𝐫^{})]=i\mathrm{}\delta ^{ij}\delta (𝐫𝐫^{}).$$
(55)
### II.2 Isolated Substrate

The vibrational modes for a semi-infinite isotropic elastic substrate, with a free surface at the $`xy`$ plane and extending to infinity in the negative $`z`$ direction, were quantized previously by Ezawa;ezawa  therefore, the details will be left out here. The displacement field can be written asfootnote1

$$𝐮_{\mathrm{sub}}(𝐫)=\underset{I}{}\sqrt{\frac{\mathrm{}}{2\rho \omega _I}}\left[b_I𝐟_I(𝐫)+b_I^{}𝐟_I^{}(𝐫)\right],$$
(56)
where $`b`$ and $`b^{}`$ are the annihilation and creation operators for the substrate phonons. The index $`I`$, like the index $`J`$ for the nanoparticle, uniquely specifies a phonon mode for the substrate. $`𝐟_I`$ are eigenfunctions of (5) subject to stress-free boundary conditions at the $`z=0`$ plane.

In what follows we will need the spectral density of the isolated substrate which is defined as

$$N_{\mathrm{sub}}(𝐫,\omega )\frac{1}{\pi }\mathrm{Im}D_{\mathrm{sub}}^{zz}(𝐫,𝐫,\omega ),$$
(57)
where $`D_{\mathrm{sub}}^{ij}(𝐫,𝐫^{},\omega )`$ is the Fourier transform of the retarded phonon Green’s function

$$D_{\mathrm{sub}}^{ij}(𝐫,𝐫^{},t)i\theta (t)[u_{\mathrm{sub}}^i(𝐫,t),u_{\mathrm{sub}}^j(𝐫^{},0)]$$
(58)
of the substrate. The spectral density at the free surface of silicon, regarding it as an isotropic elastic continuum, was calculated in Appendix B of Ref. patton and geller, . There we obtained

$$N_{\mathrm{sub}}(\omega )=C_{\mathrm{Si}}\omega ,C_{\mathrm{Si}}1.4\times 10^{46}\mathrm{cm}^2\mathrm{s}^2.$$
(59)
## III GOLDEN-RULE LIFETIMES

The relaxation rate or inverse lifetime of the perturbed eigenmodes of the nanoparticle can be calculated using Fermi’s golden rule (setting $`\mathrm{}=1`$),

$$\tau _J^1=2\pi \underset{\mathrm{f}}{}|\mathrm{f}|\delta H|\mathrm{i}|^2\delta (\omega _\mathrm{i}\omega _\mathrm{f}),$$
(60)
where the initial and final states are

$$|\mathrm{i}=a_J^{}|0\mathrm{and}|\mathrm{f}=b_I^{}|0.$$
(61)
Using Eqs. (4), (50), and (56) leads tofootnote1

$$\tau _\mathrm{J}^1=\frac{\pi K^2}{2\rho ^2}\frac{|\mathrm{\Psi }_J^z(𝐫_0)|^2}{\omega _J}\underset{I}{}\frac{|f_I^z(𝐫_0)|^2}{\omega _J}\delta (\omega _J\omega _I).$$
(62)
Noting that

$$\underset{I}{}|f_I^z(𝐫_0)|^2\delta (\omega _J\omega _I)=2\rho \omega N_{\mathrm{sub}}(\omega ),$$
(63)
we obtain (reinstating factors of $`\mathrm{}`$)

$$\tau _J^1=\frac{\pi K^2}{\mathrm{}\rho }\frac{N_{\mathrm{sub}}(\omega _J)}{\omega _J}|\mathrm{\Psi }_J^z(𝐫_0)|^2.$$
(64)
Using $`K=1.0\times 10^4`$ $`\mathrm{erg}\mathrm{cm}^2`$ and $`\rho =2.3`$ $`\mathrm{g}\mathrm{cm}^3,`$ which are appropriate (see Ref. patton and geller, ) for a weak link in Si, relaxation rates and Q factors are given in Table 1 for some low lying modes. The Q factor is defined here as $`\tau `$ divided by the period $`T`$,

$$\mathrm{Q}\frac{\tau }{T}=\frac{\mathrm{}\omega _J}{2\pi \gamma },$$
(65)
where $`\gamma \mathrm{}\tau ^1`$ is an energy width.

The values of the Q factors we obtain for the low lying modes are incredibly large, reflecting the fact that the reservoir (substrate) is extremely ineffective at absorbing energy at these low frequencies. As we will discuss below in Section VI, the lifetimes (and Q factors) for the model considered here cannot directly be compared to the experiment of Ref. yang,  without accounting for the difference in sound speeds between a solid Si substrate and a weakly bound nanoparticle cluster, as well as some other less important modifications. There we shall show that the coefficient $`C`$ in Eq. (59) should be enhanced by about 1000 before making such a comparison, which decreases the Q factors by this same factor. However, the Q factors corrected in this way are still huge, and the good agreement with the observed low-frequency DOS (see below) suggests that the Q factors of the nanoparticles studied experimentally in Ref. yang,  are also very large.

In Ref. p and g,  we used the golden rule result (64) to estimate the phonon DOS at low energies. This is achieved by replacing, in accordance with Fermi’s golden rule, each discrete mode in the isolated nanoparticle by a Lorentzian with a width given by (64). (More precisely, this amounts to approximating the energy-dependent phonon self-energy for each mode $`J`$ with its value at $`\omega =\omega _J`$, a procedure often called the quasiparticle-pole approximation.) However, this procedure is unreliable at low energies because the actual line-shapes of the broadened modes are non-Lorentzian in the tails. Nevertheless, we obtained a DOS at 3 $`\mathrm{cm}^1`$ that was only 20 times smaller than that observed.footnote3

## IV MANY-BODY THEORY OF THE DOS

### IV.1 Local DOS

To leading order in the electron-phonon interaction strength, the electronic population relaxation rate due to phonon emission (for example, as measured in Ref. yang, ) is given by Fermi’s golden rule, which states the rate (for a deformation potential electron-phonon interaction) is proportional to the square of the electron-phonon coupling strength times the phonon DOS. In a translationally invariant system the DOS does not have any position dependence, but in a nanoparticle one must distinguish between the “global” DOS (the DOS relevant for thermodynamics) and the local eigenfunction-weighted DOS, which is the one that determines phonon emission rate. We will call this position-dependent DOS the local DOS, and denote it by $`g(𝐫,\omega )`$. The precise definition of $`g(𝐫,\omega )`$ will be given below.

From a theoretical point-of-view, the quantity describing the local vibrational dynamics in the nanoparticle is the (retarded) phonon Green’s function

$$D_\mathrm{R}^{ij}(𝐫,𝐫^{},t)i\theta (t)[u^i(𝐫,t),u^j(𝐫^{},0)]_H,$$
(66)
where

$$_H\frac{\mathrm{Tr}(e^{\beta H})}{\mathrm{Tr}e^{\beta H}},$$
(67)
and with the Hamiltonian given by

$$H=H_0+\delta H.$$
(68)
Here $`H_0`$ is the Hamiltonian of the isolated nanoparticle and substrate,

$$H_0=\underset{J}{}\omega _Ja_J^{}a_J+\underset{I}{}\omega _Ib_I^{}b_I,$$
(69)
and, as mentioned in Section II, $`\delta H`$ is a harmonic spring potential given in Eq. (4).

In this section the phonon Green’s function $`D`$ always refers to the nanoparticle, and the label “nano” will be suppressed. The imaginary part of the Fourier transform of $`D_\mathrm{R}^{ij}(𝐫,𝐫,t)`$ defines the nanoparticle’s phonon spectral density

$$N^{ij}(𝐫,\omega )\frac{1}{\pi }\mathrm{Im}D_\mathrm{R}^{ij}(𝐫,𝐫,\omega ).$$
(70)
For an electron system (or any system of particles), the spectral density defined above is precisely the local DOS. However, because the elasticity equation of motion (5) is second order in time, the spectral density and DOS (both local and global) differ by a factor of $`2\rho \omega `$. In addition, the vibrational spectral density (70) is a tensor, whereas the phonon emission rate probes some coupling-constant-weighted sum of tensor elements. Because we are ascribing the observed reduction in phonon emission (in going from bulk to nanoparticle) to a reduction in the local DOS, our results are not sensitive to the precise way in which a scalar quantity is constructed from the tensor, as long as the same measure is used in both the nanoparticle and bulk. It will be most convenient to investigate the trace of the local DOS tensor. Therefore, the quantity we calculate in this paper is

$$g(𝐫,\omega )2\rho \omega \underset{i=1}{\overset{3}{}}N^{ii}(𝐫,\omega ),$$
(71)
which we shall refer to as the local DOS. $`g(𝐫,\omega )`$ characterizes the number of states per unit energy per unit volume near position $`𝐫`$. In a bulk material with Debye spectrum, (71) reduces at low frequency to

$$g(𝐫,\omega )=\frac{\omega ^2}{2\pi ^2}\left(\frac{1}{v_\mathrm{l}^3}+\frac{2}{v_\mathrm{t}^3}\right),$$
(72)
independent of $`𝐫`$. Eq. (72) is the well-known Debye formula for the vibrational DOS of a crystal.

The local DOS $`g(𝐫,\omega )`$ controls the phonon emission rate for an impurity atom sitting at position $`𝐫`$. Although the impurity locations in a real nanoparticle are assumed to be random, dopants near the surface are known to be optically inactive; hence, the experiments (including that of Ref. yang, ) do not probe the phonon DOS near the nanoparticle surface. Therefore, we introduce a particular volume-averaged DOS

$$\overline{g}(b,\omega )\frac{\underset{rb}{}d^3rg(𝐫,\omega )}{\frac{4}{3}\pi b^3},$$
(73)
which characterizes the average $`g(𝐫,\omega )`$ within a sphere of radius $`b`$. In the limit $`bR`$, in which case the local DOS is averaged over the full nanoparticle volume, we obtain the global (or thermodynamic) DOS, which, for an isolated nanoparticle, is

$$\overline{g}(R,\omega )=\frac{1}{V}\underset{J}{}\delta (\omega \omega _J).$$
(74)
Physically, we expect $`b`$ to be somewhere between $`R/2`$ and $`R`$.

### IV.2 Perturbative Calculation of the Local DOS

The retarded Green’s function (66) for the nanoparticle can be obtained by calculating the Euclidean time-ordered (or imaginary time) Green’s function defined by

$$D^{ij}(𝐫,𝐫^{},\tau )=Tu^i(𝐫,\tau )u^j(𝐫^{},0)_H.$$
(75)
In the interaction representation,

$$D^{ij}(𝐫,𝐫^{},\tau )=\frac{Tu^i(𝐫,\tau )u^j(𝐫^{},0)e^{_0^\beta \delta H(\tau ^{})𝑑\tau ^{}}_{H_0}}{e^{_0^\beta \delta H(\tau ^{})𝑑\tau ^{}}_{H_0}},$$
(76)
where the expectation values are with respect to $`H_0`$. By expanding the exponentials to leading order in the perturbation and Fourier transforming, (76) can be written as

$`D^{ij}(𝐫,𝐫^{},\omega )`$ $`=`$ $`D_0^{ij}(𝐫,𝐫^{},\omega )+{\displaystyle \underset{kl}{}}{\displaystyle D_0^{ik}(𝐫,𝐫^{\prime \prime },\omega )}`$
$`\times `$ $`\mathrm{\Pi }^{kl}(𝐫^{\prime \prime },𝐫^{\prime \prime \prime },\omega )D_0^{lj}(𝐫^{\prime \prime \prime },𝐫^{},\omega )d^3r^{\prime \prime }d^3r^{\prime \prime \prime },`$
where

$$D_0^{ij}(𝐫,𝐫^{},\omega )=\underset{J}{}\frac{\mathrm{\Psi }_J^i(𝐫)\mathrm{\Psi }_J^j(𝐫^{})}{2\rho \omega _J}\left[\frac{1}{i\omega \omega _J}\frac{1}{i\omega +\omega _J}\right]$$
(78)
is the free propagator and $`\mathrm{\Pi }^{ij}`$ is the leading order self-energy, given at zero temperature byfootnote1

$`\mathrm{\Pi }^{ij}(𝐫,𝐫^{},\omega )`$ $`=`$ $`{\displaystyle \frac{K^2}{2\rho }}{\displaystyle \underset{\mathrm{I}}{}}{\displaystyle \frac{|f_\mathrm{I}^z(𝐫_0)|^2}{\omega _\mathrm{I}}}\left[{\displaystyle \frac{1}{i\omega \omega _\mathrm{I}}}{\displaystyle \frac{1}{i\omega +\omega _\mathrm{I}}}\right]`$ (79)
$`\times `$ $`\delta ^{iz}\delta ^{jz}\delta (𝐫𝐫_0)\delta (𝐫^{}𝐫_0).`$
The $`𝐟_\mathrm{I}(𝐫)`$ are the substrate eigenfunctions discussed in Section II, and $`𝐫_0`$ is the point at which the nanoparticle is connected to the substrate. Retarded quantities are obtained by analytically continuing $`i\omega \omega +i0^+`$.

To calculate the local DOS (71) we need to solve the Dyson equation for the nanoparticle Green’s function, written symbolically as

$$D=D_0+D_0\mathrm{\Pi }D.$$
(80)
A solution to (80) can be obtained by introducing matrix representations for $`D`$, $`D_0`$, and $`\mathrm{\Pi }`$, in which case

$$D=(D_0^1\mathrm{\Pi })^1.$$
(81)
The matrix representation we use is defined by

$$O(J,J^{},\omega )\underset{ij}{}_Vd^3rd^3r^{}\mathrm{\Psi }_J^i(𝐫)\mathrm{\Psi }_J^{}^j(𝐫^{})O^{ij}(𝐫,𝐫^{},\omega ),$$
(82)
where $`O=D`$, $`D_0`$, or $`\mathrm{\Pi }`$. In (82) the integration is over the volume $`V`$ of the nanoparticle, and the $`𝚿_J(𝐫)`$ are the nanoparticle eigenfunctions. The inverse transformation is

$$O^{ij}(𝐫,𝐫^{},\omega )=\underset{JJ^{}}{}\mathrm{\Psi }_J^i(𝐫)O(J,J^{},\omega )\mathrm{\Psi }_J^{}^j(𝐫^{}).$$
(83)
A nanoparticle with a diameter of 10 nm has approximately 8,000 atoms in it, thus, there are roughly 24,000 acoustic vibrational modes. By knowing the number of modes, a Debye energy can be defined: the Debye energy is the energy at which there are 24,000 elasticity-theory modes that lie below in energy. For our nanoparticle, the Debye energy is about 240 cm<sup>-1</sup>. The Debye energy cutoff truncates the Hilbert space, which leads to finite-size matrices.footnote4  This enables every mode $`J`$ of the nanoparticle to be included in the calculation of the Green’s function (81).

## V RESULTS

In this section we present our results for the phonon DOS in a 10 nm Si nanoparticle, obtained by solving the Dyson equation (80) for the phonon Green’s function, as explained above. As we have discussed, the DOS $`g(𝐫,\omega )`$, defined in (71), is a local quantity that varies with position within the nanoparticle, and, as mentioned in the previous section, the quantity we are interested in is $`\overline{g}(b,\omega )`$, which is $`g(𝐫,\omega )`$ averaged over a sphere of radius $`b`$ centered at the center of the nanoparticle. Because we have found no significant dependence of $`\overline{g}(b,\omega )`$ on $`b`$, for the physically relevant values of $`b`$ ($`R/2<b<R`$), we plot $`\overline{g}(R,\omega )`$. As stated above in Section IV.1, $`\overline{g}(R,\omega )`$ is the global phonon DOS in the nanoparticle.

For simplicity we assume both the nanoparticle and the substrate to be made of Si; this allows us to use the surface spectral density (59) calculated in Appendix B of Ref. patton and geller, , where Si is treated as an isotropic elastic continuum with longitudinal and transverse sound velocities

$`v_\mathrm{l}=8.5\times 10^5\mathrm{cm}\mathrm{s}^1,`$
$`v_\mathrm{t}=5.9\times 10^5\mathrm{cm}\mathrm{s}^1,`$ (84)
and mass density $`\rho =2.3\mathrm{g}\mathrm{cm}^3`$. In the final section of this paper, where we compare our results to the experiment of Ref. yang, , we will introduce an important correction to account for the differences between a solid Si substrate and a nanoparticle cluster.

In Fig. 2 the global DOS $`\overline{g}(R,\omega )`$ of a 10 nm diameter nanoparticle is given up to 100 cm<sup>-1</sup>. The modes above 100 cm<sup>-1</sup> were included in the calculation, but the long-wavelength approximation of elasticity theory becomes invalid at high energy. Thus, only the lower part of the spectrum is shown. Fig. 3 shows the low-energy phonon DOS up to about 15 cm<sup>-1</sup>. The large peak on the right side is the well-known Lamb mode. The phonon DOS at 3 cm<sup>-1</sup> is approximately $`4.9\times 10^{10}`$ states per wavenumber per cm<sup>3</sup>.

## VI COMPARISON WITH EXPERIMENT

In this section we compare our results with the experiment of Ref. yang, , where the one-phonon emission rate (and, therefore, the phonon DOS at 3 cm<sup>-1</sup>) in a cluster of Y<sub>2</sub>O<sub>3</sub> nanoparticles was observed to be $`8.2\times 10^3`$ times that in bulk Y<sub>2</sub>O<sub>3</sub>. In particular, the excited $`{}_{}{}^{5}\mathrm{D}_{1}^{}(\mathrm{II})`$ state of Eu<sup>3+</sup> in the nanoparticles had a phonon-emission lifetime of 27$`\mu `$s, compared with a bulk value of 221 ns. In order to make a comparison of our results to that of the experiment, two modifications of our calculation have to be performed.

In our model, the cluster of nanoparticles has been replaced by a solid substrate. However, the spectral density (59) of the substrate, which at long wavelengths, is determined by the sound speeds and mass density of Si, is much smaller than that of the nanoparticle cluster. Treating the long-wavelength modes of the cluster with elasticity theory (or, even simpler, approximating the random cluster by an ordered cubic lattice), shows that the spectral density (59) should be replaced by (the subscript “cl” referring to cluster)

$$N_{\mathrm{cl}}(\omega )=C_{\mathrm{cl}}\omega ,$$
(85)
wherefootnote5

$$C_{\mathrm{cl}}=\frac{v_{\mathrm{Si}}^3\rho _{\mathrm{Si}}}{v_{\mathrm{cl}}^3\rho _{\mathrm{cl}}}C_{\mathrm{Si}}.$$
(86)
Here $`v_{\mathrm{cl}}`$ is a characteristic sound speed in the cluster and $`\rho _{\mathrm{cl}}`$ is its mass density. The $`1/v_{\mathrm{cl}}^3`$ dependence in (86) comes from the well-known velocity dependence of the Debye DOS, and the $`1/\rho _{\mathrm{cl}}`$ factor comes from the definition of spectral density \[see discussion following Eq. (70)\]. Approximating the cluster by an ordered cubic array with lattice constant $`d`$ (the nanoparticle diameter), yields

$$v_{\mathrm{cl}}\sqrt{\frac{K}{M}}d$$
(87)
and

$$\rho _{\mathrm{cl}}\frac{\pi }{6}\rho _{\mathrm{Si}},$$
(88)
where $`K`$ is the effective spring constant connecting the nanoparticles, given after Eq. (64), and

$$M=\frac{4}{3}\pi \rho _{\mathrm{Si}}\left(\frac{d}{2}\right)^3$$
(89)
is the mass of one nanoparticle. Using $`d=13`$ nm (the mean diameter in Ref. yang, ), we obtain an enhancement factor

$$\frac{C_{\mathrm{cl}}}{C_{\mathrm{Si}}}=9.8\times 10^2.$$
(90)
This factor increases the nanoparticle DOS (at frequencies below the Lamb mode) by nearly three orders of magnitude.

There are several other marginally important corrections, most of which will be ignored, and one that will be included just for completeness: In the model analyzed above, the nanoparticle was connected to its surroundings by only a single contact point, whereas a nanoparticle in a cluster most likely has more than one connection. As the number of contacts increases, this simply scales the DOS (away from the peaks) linearly with the number of contact points. Conservatively, we expect that the multiple contact points present in the real system will increase the decay rate of the nanoparticles’ modes, and hence the phonon DOS well below the Lamb mode, by a factor of two.

The following additional corrections have also been considered and were found not to be significant (and are not included in our final results):

(i) The actual experiment of Ref. yang,  was done on an ensemble of nanoparticles with mean diameter 13 nm and standard deviation of 5 nm. To understand the effects of this size distribution, we have calculated the DOS at 3 cm<sup>-1</sup> averaged over a Gaussian distribution of diameters centered at 10 nm. Even for very wide distributions (standard deviation up to 8 nm), the ensemble averaged DOS at 3 cm<sup>-1</sup> is increased by no more than a factor of two.footnote2  In addition, the correction for re-centering the size distribution from 10 nm to the experimentally observed 13 nm average size, leads to corrections of order unity.

(ii) Our calculations were done for a nanoparticle and substrate made of Si, while the experiment was done on Y<sub>2</sub>O<sub>3</sub> nanoparticles, which, of course, has different mass density and sound velocities. The differences in sound speeds and mass only shift around the modes of the nanoparticle by a small amount; this has the same effect as a small change in the diameter of the nanoparticles, which we have found to be negligible. As for the substrate (or more precisely the replacement of the substrate with a cubic lattice of nanoparticles), the change in mass density does effect the spectral function (59) by changing the velocity (87) and mass density (88), but this change is only of order unity.

(iii) The nanoparticles of Ref. yang,  were immersed in He, either liquid (for T$`<`$4.21 K) or gas (T$`>`$4.21 K). However, the results were found not to change through the liquid-gas transition, presumably because of the large sound-speed mismatch between superfluid He and Y<sub>2</sub>O<sub>3</sub>. Therefore, we have ignored the presence of He in our theory.

(iv) The experiment of Ref. yang,  was done at temperatures between 1.5 and 10 K (excluding the interval 2.17 to 4.21 K), whereas our calculations assume zero temperature. The effect of finite temperature is to stimulate phonon emission into the bath (i.e. substrate or cluster). However, this is not important until the Bose distribution function of the bath at the Lamb mode frequency (approximately 13 cm<sup>-1</sup>) becomes of order unity, which does not occur until the temperature exceeds about 18 K.

Taking into account the first two modifications, and ignoring the rest, we obtain a 3 cm<sup>-1</sup> DOS given by

$`\overline{g}(R,3\mathrm{c}\mathrm{m}^1)`$ $`=`$ $`4.9\times 10^{10}{\displaystyle \frac{\mathrm{states}}{\mathrm{wavenumber}\mathrm{cm}^3}}`$ (91)
$`\times `$ $`980\times 2`$
$`=`$ $`9.6\times 10^{13}{\displaystyle \frac{\mathrm{states}}{\mathrm{wavenumber}\mathrm{cm}^3}}.`$
The DOS of bulk Si is given by

$$\underset{\lambda }{}\frac{ϵ^2}{2\pi ^2\mathrm{}^3v_\lambda ^3}=3.9\times 10^{15}E^2\frac{\mathrm{states}}{\mathrm{wavenumber}\mathrm{cm}^3},$$
(92)
where $`E`$ is the energy in wavenumbers. Thus, the theoretical ratio of the DOS of the nanoparticle to that of the bulk material at 3 cm<sup>-1</sup> is $`2.7\times 10^3`$. As stated above, the experimental ratio of nanoparticle to bulk DOS was found to be approximately $`8.2\times 10^3`$. The agreement between our theory and the experiment of Ref. yang,  is excellent considering the simplicity and robustness of our model. We conclude that the low-energy phonon DOS observed in Ref. yang,  is consistent with our enviromental broadening mechanism.

## VII CONCLUSIONS

Motivated by an interesting experimentyang  measuring the low-energy phonon DOS in a insulating nanoparticle, we have thoroughly investigated a simplified model of a single nanoparticle weakly coupled to its environment, a semi-infinite substrate. The environmental interactions were found to significantly affect the DOS at energies below the Lamb mode.

Additionally, we have used the results of our model to predict the effect of environmental interaction in a cluster of nanoparticles like that studied in Ref. yang, . Although it is necessary to estimate the value of several quantities appearing in the model, we believe that we can do this accurately enough to obtain a final result that is correct at the order-of-magnitude level, with no free parameters. Because our results for the 3 cm<sup>-1</sup> DOS is only about three times smaller than that observed in Ref. yang, , our broadening mechanism and the resulting phonon spectrum is clearly consistent with that experiment.

## VIII ACKNOWLEDGMENTS

This work was supported by the National Science Foundation under CAREER Grant No. DMR-0093217, and by the Research Corporation. Is is a pleasure to thank Steve Lewis, Vadim Markel, Richard Meltzer, Mark Stockman, and Ho-Soon Yang for useful discussions, and Patrick Sprinkle for help with the numerics. We would especially like to thank Bill Dennis for his help with every aspect of this work.


format: text

text:
# Yang-Lee zeros for a nonequilibrium phase transition

## 1 Introduction

The investigation of nonequilibrium systems is of great importance since most phenomena in nature take place under nonequilibrium conditions. In this field, as in equilibrium physics, phase transitions are particularly interesting. However, when dealing with nonequilibrium systems one cannot utilize such a well-established theoretical framework as in equilibrium statistical mechanics. Therefore it is interesting to investigate which concepts of equilibrium physics can be transfered to nonequilibrium systems.

An important concept of equilibrium statistical physics is the theory of Yang and Lee  for the emergence of nonanalytic behaviour at phase transtitions. Only recently these ideas have been applied to an integrable nonequilibrium model . In the present work we use similar techniques to investigate directed percolation (DP)  as a paradigm for nonintegrable systems far from equilibrium. We show that it is possible to transfer the ideas of Yang and Lee to DP and that this method does indeed provide information about universal properties of the phase transition.

### 1.1 Yang and Lee’s theory in equilibrium statistical physics

We shall shortly sketch the ideas of complex zeros of the partition function in equilibrium statistical mechanics. Taking the Ising model without an external magnetic field as an example, we consider a system of $`N`$ spins $`s_i`$ ($`s_i=\pm 1,i=1,\mathrm{},N`$) on a lattice where the spins interact with their nearest neighbours, denoted by $`<i,j>`$. The total energy $`(\left\{s_i\right\})`$ of a configuration $`\left\{s_i\right\}`$ of the spins is given by $`(\left\{s_i\right\})=J_{<i,j>}s_is_j`$. The system is said to be in thermal equilibrium at temperature $`T`$ when the probability distribution $`P(\left\{s_i\right\})`$ to find the system in configuration $`\left\{s_i\right\}`$ is given by the stationary Gibbs ensemble $`P(\left\{s_i\right\})\mathrm{}^{\beta (\left\{s_i\right\})},\beta =k_\mathrm{B}T`$, where $`k_\mathrm{B}`$ is Boltzmann’s constant. In that case the statistical properties of the system are fully determined by the canonical partition function

$`Z_N(T)`$ $`=`$ $`{\displaystyle \underset{\left\{s_i\right\}}{}}\mathrm{}^{\beta (\left\{s_i\right\})}={\displaystyle \underset{\left\{s_i\right\}}{}}x^{r(\left\{s_i\right\})}\mathrm{with}`$ (2)
$`x=\mathrm{}^{\beta J}>0,r(\left\{s_i\right\})={\displaystyle \underset{<i,j>}{}}s_is_j.`$
For finite system size $`Z_N`$ can be expressed as the ratio of two polynomials in $`x`$ with integer coefficients. The free energy per lattice site is obtained from $`Z_N`$ through

$`f=\underset{N\mathrm{}}{lim}\left({\displaystyle \frac{1}{\beta N}}\mathrm{ln}Z_N\right)`$ (3)
The hallmark of phase transitions is the appearance of nonanalyticities. Since in Equation (3) singularities occur as zeros of the partition function $`Z_N`$, these are good candidates for indicating a phase transition. However, $`Z_N`$ has no zeros for a positive real temperature so that there is no phase transition for finite $`N`$. This led Yang and Lee to study the complex zeros of $`Z_N`$555For a field-driven phase transition these zeros are denoted as Yang-Lee zeros while in the case of a temperature-driven phase transition the term Fisher zeros is used., which in this case correspond to complex temperatures, and examine their behaviour as the system size grows. They argued that in the thermodynamic limit $`N\mathrm{}`$ in (3), if the model exhibits a phase transition, an ever increasing number of complex zeros accumulate in the vicinity of a point on the positive real axis, which corresponds to a physical temperature. The zeros thereby induce singularities at the critical point in the thermodynamic limit, thus explaining the crossover to nonanalyticity at the transition in the limit of an infinite system size. Such a behaviour has been observed for many equilibrium models. Moreover, several features of the distribution of zeros have been related to universal properties of the system under consideration .

### 1.2 Directed percolation, a process away from equilibrium

Since nonequilibrium systems do not obey the stationary Gibbs ensemble they cannot be handled using the partition function of equilibrium statistical physics. Instead they are described by the time-dependent probability distribution $`P_t(\left\{s_i\right\})`$ which has to be derived from the master equation. In most cases this cannot be done exactly. Nevertheless, the concept of universality, well-known from equilibrium statistical physics, proves to be suitable for nonequilibrium phase transitions as well. In this context models exhibiting a transition from a fluctuating active phase to a nonfluctuating inactive phase (absorbing state) have been extensively studied . These models are used to describe spreading processes as e.g. forest fires, where a spreading agent can either spread over the entire system or die out after some time. The two phases are seperated by a nonequilibrium phase transition. The most prominent universality class of transitions into an absorbing state is that of directed percolation. Directed percolation is an anisotropic variant of ordinary percolation in which activity can only percolate along a given direction in space. Regarding this direction as a temporal degree of freedom, DP can be interpreted as a dynamical process. Directed percolation emerges in a variety of physical problems ranging from catalytic reactions on surfaces  to spatio-temporal intermittency in magnetic fluids .

### 1.3 Realizations of DP

Simple realizations of DP are directed bond (DPb) and directed site percolation (DPs) on a tilted square lattice (see Figure 1). In directed bond percolation the bonds are conducting with probability $`p`$ and non-conducting with probability $`1p`$. In this model sites at time $`t>0`$ are activated by directed paths of conducting bonds, originating from active sites at time $`t=0`$. A cluster consists of all sites that are connected by such paths of conducting bonds to active sites in the initial state. In directed site percolation on the other hand all bonds are conducting while the sites themselves can be either permeable ($`p`$) or blocked ($`1p`$). Activity can spread from permeable site to permeable site. A cluster is formed by permeable sites that are connected to active sites at time $`t=0`$ by a directed path of bonds that only connect permeable sites.

The order parameter which characterizes the phase transition is the probability $`P(\mathrm{})`$ that a randomly chosen site belongs to an infinite cluster. For $`p>p_\mathrm{c}`$ this probability is finite whereas it vanishes for $`pp_\mathrm{c}`$. Close to the phase transition $`P(\mathrm{})`$ is known to vanish algebraically as $`P(\mathrm{})(pp_\mathrm{c})^\beta `$, where $`\beta `$ is an universal critical exponent. In addition, the DP process is characterized by a spatial correlation length $`\xi _{}`$ (perpendicular to time) and a temporal correlation length $`\xi _{||}`$. As $`p`$ approaches $`p_\mathrm{c}`$ these length scales are known to diverge as

$`\xi _{}|pp_\mathrm{c}|^\nu _{},\xi _{||}|pp_\mathrm{c}|^{\nu _{||}}`$ (4)
with the critical exponents $`\nu _{}`$ and $`\nu _{||}`$. The scaling behaviour (4) implies that DP is invariant under the scaling transformations

$`x\mathrm{\Lambda }x,t\mathrm{\Lambda }^zt,(pp_\mathrm{c})\mathrm{\Lambda }^{1/\nu _{}}(pp_\mathrm{c}),`$ (5)
where $`x`$ is the position and $`z=\nu _{||}/\nu _{}`$ is the so-called dynamic exponent. Numerical estimates for the values of $`p_\mathrm{c}`$ and the critical exponents in 1+1 dimensions are $`p_\mathrm{c}=0.6447001(1)`$, $`\beta =0.27649(4)`$, $`\nu _{}=1.096844(14)`$ and $`\nu _{||}=1.733825(15)`$  (bond percolation) and $`p_\mathrm{c}=0.705489(4)`$  (site percolation).

Although DP can be defined and easily simulated, it is one of the very few systems for which – even in one spatial dimension – no analytical solution is known, suggesting that DP is a non-integrable process. In fact, the values of the percolation threshold and the critical exponents are not simple numerical fractions but seem to be irrational instead.

### 1.4 Yang and Lee zeros for directed percolation

To apply the idea of Yang-Lee zeros to DP, we consider the order parameter in a finite system as a function of the percolation probability $`p`$ in the complex plane. This can be done by studying the time-dependent survival probability $`P(t)`$, which is defined as the probability that a cluster generated in a single site at time $`t=0`$ survives up to time $`t`$ (or even longer). For a finite system $`P(t)`$ can be expressed as a polynomial in $`p`$ (see below). Note that $`lim_t\mathrm{}P(t)`$ and the order parameter $`P(\mathrm{})`$ coincide.

The survival probability and the partition function (2) show a similar behaviour in many respects. For finite systems $`P(t)`$ and $`Z_N`$ do not have relevant zeros in the physical region of the control parameter $`0p1`$ and $`0<T<\mathrm{}`$ although the phase transition is marked by a vanishing $`P(t)`$ and $`Z_N`$ at the critical points in the limit of infinite systems. At those points both functions exhibit nonanalytic behaviour. The zeros of $`P(t)`$ and $`Z_N`$ are generated by polynomials with real integer coefficients, thus the zeros come in complex-conjugate pairs.

According to Yang and Lee the increasing system size is accompanied by an approach of some of the complex zeros of $`Z_N`$ to the physical region of the control parameter, while the accumulation point marks the critical point. We show that the complex zeros of $`P(t)`$ approach the critical value $`p_\mathrm{c}`$ on ’trajectories’ for increasing time $`t`$ and that the distance between the zeros and the critical point is related to the critical exponent for the temporal correlation length.

### 1.5 Determining the survival probability

In directed bond (site) percolation the survival probability $`P(t)`$ is given by the sum over the weights of all possible configurations of bonds (sites) for which the process survives at least up to time $`t`$. Each conducting bond (permeable site) contributes to the weight with a factor $`p`$, while each non-conducting bond (blocked site) contributes with a factor $`1p`$. However, the states of those bonds (sites) which do not touch the actual cluster are irrelevant as they do not contribute to the survival of the cluster. Therefore, it is sufficient to consider the sum over all possible clusters $`𝒞`$ of bonds (sites) connected to the origin. Each cluster is weighted by the contributions of the conducting bonds (permeable sites) belonging to the cluster and the non-conducting bonds (blocked sites) belonging to its hull. Roughly speaking, the hull surrounds the cluster. More precisely the hull of a cluster is the set of non-conducting bonds (blocked sites) that would contribute to the cluster if they were conducting (permeable) (see Figure 2). Thus the survival probability can be expressed as

$$P(t)=\underset{𝒞}{}p^n(1p)^m,$$
(6)
where the sum runs over all clusters reaching the horizontal row at time $`t`$. For each cluster $`n`$ denotes the number of its bonds (sites), while $`m`$ is the number of bonds (sites) belonging to its hull. Note that in this sense the hull does not include bonds connecting sites at time $`t`$ and $`t+1`$ (DPb) or sites at time $`t+1`$ (DPs) since the cluster may survive even longer. Summing up all weights in Equation (6) one obtains a polynomial. As can be verified, the first few polynomials for the survival probability in the directed bond percolation process are

$`P(0)=1`$ (7)
$`P(1)=2pp^2`$
$`P(2)=4p^22p^34p^4+4p^5p^6`$
$`P(3)=8p^34p^410p^53p^6+18p^7+5p^830p^9+24p^{10}8p^{11}+p^{12}`$
$`P(4)=16p^48p^524p^68p^7+6p^8+84p^929p^{10}62p^{11}120p^{12}`$
$`+244p^{13}+75p^{14}470p^{15}+495p^{16}268p^{17}+83p^{18}14p^{19}+p^{20}.`$
As $`t`$ increases the number of cluster configurations grows rapidly, leading to complicated polynomials with very large coefficients (e.g. for $`t=15`$ the largest coefficient for bond percolation is of order $`10^{44}`$).

The paper is structured as follows. In Section 2 we investigate if the distribution of zeros is totally dependent on microscopic details of the underlying model or whether it also contains universal features. We ask whether the distribution is related to some of the critical exponents which characterize the phase transition of directed percolation. Section 3 presents certain non-trivial points where the polynomials $`P(t)`$ for bond percolation can be calculated exactly for all times. Section 4 ends with a conclusion. In A we show that the value of the survival probability can be computed by Monte Carlo simulations even if $`p`$ is a complex number. The first coefficients of the polynomials $`P(t)`$ for bond percolation are discussed in B.

## 2 Universal properties

The distribution of the zeros of $`P(t)`$ (from $`t=8`$ to $`t=15`$) in the neighbourhood of the critical point is shown in Figure 3.

Away from the critical point the appearance of the distributions for DPb and DPs is quite different. This implies that in this region the distribution of zeros depends strongly on microscopic details and hence is non-universal. However, there is also a general feature which can be observed in both cases. The innermost zeros approach the critical point on ‘trajectories’ as $`t`$ increases (see Figure 4).

We applied a standard Bulirsch-Stoer (BST) acceleration algorithm  to the set of zeros of each enumerated ‘trajectory’ of Figure 4 to determine $`lim_t\mathrm{}\mathrm{Re}\left(p^{\mathrm{zero}}(t)\right)`$ and $`lim_t\mathrm{}\mathrm{Im}\left(p^{\mathrm{zero}}(t)\right)`$. The results are listed in Table 1.

Although we calculated the zeros only for small systems (until $`t=15`$) the extrapolants accord fairly well with the numerical values of the percolation threshold666The convergence for site percolation is slower than for bond percolation since the order of the polynomials $`P(t)`$ for DPs is smaller than for DPb. $`p_\mathrm{c}=0.6447001(1)`$ (bond percolation) and $`p_\mathrm{c}=0.705489(4)`$ (site percolation). Thus a similar scenario as observed by Yang and Lee for equilibrium phase transitions proves to be suitable for the nonequilibrium phase transition of DP. As time increases the zeros of the survival probability $`P(t)`$ approach the real axis between $`p=0`$ and $`p=1`$. The accumulation point is the critical point.

So far we have shown that the zeros of $`P(t)`$ provide information about the existence of the phase transition and about the critical value $`p_\mathrm{c}`$. However, the value of $`p_\mathrm{c}`$ is non-universal and depends on the particular realization of DP. On the contrary universal features are independent of the microscopic details of the underlying process. On each ‘trajectory’ of Figure 4 we calculated the distance between the zeros and the critical point $`d(t)=|p^{\mathrm{zero}}(t)p_c|`$ with the values of $`p_\mathrm{c}`$ as given in section 1.3. According to Equation (4) and simple scaling arguments we expect $`d(t)`$ to decrease as $`d(t)t^{1/\nu _{||}}`$. Application of the BST algorithm yields results which support this claim, as shown in Table 2.

The extrapolants are in fairly good agreement with the numerical value of $`1/\nu _{||}=0.57676`$, even for small systems. This means that universal properties of the phase transition of directed percolation are encoded in the complex zeros of the survival probability $`P(t)`$. Figure 5 shows the used data and the asymptotic power law $`d(t)t^{1/\nu _{||}}`$ with the numerical value of $`\nu _{||}`$.

## 3 Exact results

In this section we will address a particularly surprising observation, namely the existence of certain points on the real axis where the polynomials $`P(t)`$ for bond percolation can be solved exactly for all values of $`t`$. Beside the trivial points $`p=0`$ (where $`P(t)=\delta _{t,0}`$) and $`p=1`$ (where $`P(t)=1`$) we find a $`t`$-independent zero at $`p=2`$ and, even more surprisingly, a very simple solution if $`p`$ is equal to one of the Golden Ratios $`(1\pm \sqrt{5})/2`$. The Golden Ratios are the roots of the quadratic equation $`p^2=p+1`$ and play an important role not only in number theory  but also in other fields ranging from chaotic systems  to arts . Although these special points are located outside the physically accessible region $`0p1`$, their existence may be helpful for further investigations of the polynomials $`P(t)`$.

### 3.1 Time-independent zero at $`p=2`$

For $`p=2`$ and $`t1`$ all polynomials $`P(t)`$ vanish identically. This can be shown as follows. Let us consider the probability $`R(t)`$ that a cluster dies out at time $`t`$, i.e., the row at time $`t`$ is the last row reached by a cluster. Obviously $`R(t)`$ is related to the survival probability by

$$R(t)=P(t)P(t+1).$$
(8)
Clearly, $`R(t)`$ can be expressed as a weighted sum over the same set of clusters as in (6). However, in the present case the weights differ from those in Equation (6) by the number of non-conducting bonds in the clusters hull between $`t`$ and $`t+1`$ since it is now required that all sites at time t+1 are inactive (see Figure 6). This means that $`R(t)`$ can be expressed as

$$R(t)=\underset{𝒞}{}p^n(1p)^m(1p)^{2k},$$
(9)
where $`n`$, $`m`$ and $`𝒞`$ have the same meaning as in (6) and $`k`$ is the number of active sites in the horizontal row at time $`t`$. Obviously, for $`p=2`$ the additional factor $`(1p)^{2k}`$ drops out so that $`P(t)=R(t)`$ for all values of $`t`$. Moreover, $`R(0)=P(0)=1`$ (for $`p=2`$). Combining these results with Equation (8) we arrive at $`P(t)=0`$ for $`t>0`$, which completes the proof.

### 3.2 Exact solution for $`p`$ at the Golden Ratio

For $`p=(1\pm \sqrt{5})/2`$ we find that the survival probability ‘oscillates’ between two different values, namely

$$P(t)=\{\begin{array}{cc}1\hfill & \text{if }t\text{ is even}\hfill \\ \frac{\pm \sqrt{5}1}{2}\hfill & \text{if }t\text{ is odd.}\hfill \end{array}$$
(10)
To prove this result, we first verify that (10) is indeed satisfied for $`t=0`$ and $`t=1`$. Then we show that

$$P(t)=P(t2)\text{for }t2\text{and}p=(1\pm \sqrt{5})/2.$$
(11)
However, instead of analyzing the survival probability directly, it turns out to be more convenient to consider the complementary probability $`Q(t)=1P(t)`$ that a cluster does not survive until time $`t`$. Obviously, $`Q(t)`$ is the sum over the weights of all clusters which do not reach the horizontal row at time $`t`$, i.e., we impose the boundary condition $`s_{0,t}=s_{1,t}=\mathrm{}=s_{t,t}=0`$. Depending on the states of the two sites $`s_{0,t1}`$ and $`s_{0,t2}`$ at the left edge of the clusters, this set of clusters may be separated into three different subsets, namely,

(a) a subset where $`s_{0,t1}=s_{0,t2}=1`$,

(b) a subset where $`s_{0,t1}=0`$ and $`s_{0,t2}=1`$, and

(c) a subset where $`s_{0,t1}=s_{0,t2}=0`$.

Next we show that the weights of the clusters in the subsets (a) and (b) cancel each other. To this end we note that the weighted sum $`\widehat{Q}(t)`$ over all clusters in subsets $`(a)`$ and $`(b)`$ may be decomposed into two independent factors $`\widehat{Q}(t)=\widehat{Q}_1\widehat{Q}_2`$, where $`\widehat{Q}_1`$ depends only on the state of the three bonds between the sites $`s_{0,t2}`$, $`s_{0,t1}`$, $`s_{0,t}`$, and $`s_{1,t}`$ (inside the box in Figure 7), while $`\widehat{Q}_2`$ accounts for all other relevant bonds. Obviously, the first factor is given by

$$\widehat{Q}_1^{(a)}=p(1p)^2,\widehat{Q}_1^{(b)}=1p,$$
(12)
while $`\widehat{Q}_2`$ takes the same value in both subsets. Thus, if $`p`$ is given by the Golden Ratio, we obtain $`\widehat{Q}_1^{(a)}+\widehat{Q}_1^{(b)}=\widehat{Q}_1(t)=0`$ and therefore the weights of subsets (a) and (b) cancel each other. Consequently, all remaining contributions to $`Q(t)`$ come from the clusters in subset (c) where the sites $`s_{0,t1}`$ and $`s_{0,t2}`$ are inactive. Now we can iterate this procedure by successively considering the sites $`s_{j,t1}`$ and $`s_{j,t2}`$ from the left to the right, where $`j=1\mathrm{}t2`$. In this way it can be shown that all these sites have to be inactive as well. Therefore, the only surviving contributions are those in which the entire row of sites at $`t2`$ is inactive, implying that $`Q(t)=Q(t2)`$. The proof of Equation (10) then follows by induction.

## 4 Conclusions

In this paper we have investigated the applicability of Yang-Lee theory to the nonequilibrium phase transition of directed percolation. In Section 1 we developed an idea of how to transfer the concepts of Yang and Lee to the phase transition of DP by studying the complex zeros of the survival probability. In Section 2 we showed that a similar scenario as observed for equilibrium phase transitions is also suitable for DP, namely complex zeros approaching the critical point for increasing system size. Moreover, we could extract the value of the critical exponent for the temporal correlation length from the distribution of zeros. Hence the zeros encode univerals properties of the phase transition of DP. Section 3 presented exact results on the survival probability of bond percolation which may be helpful for further investigations. More precisely, we proved that there are certain non-trivial values of $`p`$ for which the values of the polynomials $`P(t)`$ can be calculated exactly for all times.

This work was supported by the DAAD/CAPES within the German-Brazilian cooperation project PROBRAL – “Rigorous Results in Nonequilibrium Statistical Mechanics and Nonlinear Physics”.

## Appendix A MC-Simulations

Here we demonstrate that the value of the survival probability for a complex percolation probability $`p`$ as well as the polynomials $`P(t)`$ are accessible by computer simulations.

Let us consider the definition of the survival probability as a sum over configurations of surviving clusters in Equation (6). Formally this expression can be rewritten as

$$P(t)=\underset{𝒞}{}p^n(1p)^m=\underset{𝒞}{}q^n(1q)^m\underset{f_{n,m}}{\underset{}{\frac{p^n(1p)^m}{q^n(1q)^m}}}$$
(13)
where $`n`$ and $`m`$ denote the number of bonds of the cluster and of its hull respectively. Therefore, instead of simulating the system at a given $`p`$, we may simulate it using a different percolation probability $`q`$ reweighting each cluster by the factor $`f_{n,m}`$. While $`q(0,1)`$ still has to be a real number, $`p`$ is no longer restricted to be real, it can be any complex number. Using this reweighting technique it is in principle possible to access the entire complex plane by numerical simulations. For any finite $`t`$ such a simulation is stable and will converge to the correct result. But even for small $`t`$ the convergence time can be very long, limiting the range of applications.

Using the reweighting technique it is also possible to approximate the coefficients of the polynomial $`P(t)=_ka_kp^k`$ by considering $`p`$ as a free parameter and expanding the term $`(1p)^m`$ in Equation (13). The coefficients are then given by

$$a_k=\underset{𝒞}{}q^n(1q)^m\frac{C(m,kn)(1)^{kn}}{q^n(1q)^m},$$
(14)
where $`q(0,1)`$ is again a free parameter and

$$C(m,kn)=\{\begin{array}{cc}\frac{m!}{(kn)!(mkn)!}& \text{ if }0knm\\ 0& \text{ otherwise }.\end{array}$$
(15)
However, in most cases the direct construction of the polynomials using symbolic algebra turns out to be more efficient.

## Appendix B First coefficients of the polynomials $`P(t)`$ in the limit $`t\mathrm{}`$

As can be seen in Equation (7) the first non-vanishing coefficient of the polynomial for the survival probability $`P(t)=_na_np^n`$ for DPb is always a power of $`2`$. This contribution corresponds to the configurational weight of a surviving path without branches. The purpose of this appendix is to point out that even the following coefficients can be computed exactly, provived that $`t`$ is large enough. More specifically, we conjecture that

$$a_{t+k}=2^{tk1}q_k(t),$$
(16)
where

$$q_k(t)=\underset{m=0}{\overset{M}{}}b_m^{(k)}t^m,M=\{\begin{array}{cc}k/2\hfill & \text{ if }k\text{ even}\hfill \\ (k1)/2\hfill & \text{ if }k\text{ odd}\hfill \end{array}$$
(17)
is a finite polynomial in $`t`$. Using the explicit expressions for $`P(t)`$ up to $`t=15`$ we find that the first eight polynomials read

$`q_0\left(t\right)`$ $`=`$ $`2`$
$`q_1\left(t\right)`$ $`=`$ $`2`$
$`q_2\left(t\right)`$ $`=`$ $`42t`$
$`q_3\left(t\right)`$ $`=`$ $`2t`$
$`q_4\left(t\right)`$ $`=`$ $`16+3t+t^2`$ (18)
$`q_5\left(t\right)`$ $`=`$ $`32+t+3t^2`$
$`q_6\left(t\right)`$ $`=`$ $`232+{\displaystyle \frac{64}{3}}t+3t^2{\displaystyle \frac{1}{3}}t^3`$
$`q_7\left(t\right)`$ $`=`$ $`808{\displaystyle \frac{34}{3}}+5t^2{\displaystyle \frac{5}{3}}t^3`$
We conjecture that the leading coefficient of these polynomials is given by

$$b_M^{(k)}=\{\begin{array}{cc}2(1)^M/M!\hfill & \text{ if }k\text{ even }\hfill \\ 2(k1)(1)^M/M!\hfill & \text{ if }k\text{ odd }.\hfill \end{array}$$
(19)
This implies that the first non-vanishing coefficients of the polynomial $`P(t)`$ grow in such a way that their limit $`lim_t\mathrm{}\frac{a_{t+k}}{2^tt^M}`$ for fixed $`k`$ is well-defined. Summing up these contributions we obtain for the survival probability

$`P(t)`$ $``$ $`(2p)^t{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(p^2t/4)^j}{j!}}\left(1+p(j1/2)\right)`$
$`=`$ $`{\displaystyle \frac{1}{4}}(2p)^te^{p^2t/4}(42pp^3t).`$
Physically this expression corresponds to loop-free graphs and serves as an approximation close to $`p=0`$. The convergence radius of Equation (B) determined from $`2p=e^{p^2/4}`$ is $`|p|<0.53744`$. Note that none of the non-trivial roots computed up to $`t=15`$ lies inside this radius. We conjecture that this might be true for any $`t`$.

## References