text:
# Generation of entangled 3D localized quantum wave-packets via optical parametric amplification

## I Introduction

The quantum description of diffracting non-monochromatic beams has been the subject of intense research (see e.g. references Chiao91 ; Huttner92 ; Kennedy88 ; Drummond99 ; Dung98 ). Various approaches have been proposed in order to provide a consistent quantum picture of nonlinear optical processes and related phenomena, like entanglement and photon localization. Several important issues have been discussed in the literature, as equal-time or equal-space commutation rules, or paraxial approximations, and different “schemes” have been formulated. Kozlov01 ; Haus90 ; Deutsch91 ; Crosignani95 ; Hagelstein96  Most of them are based on the second quantization of the relevant classical equations, or on a three dimensional plane wave expansion, however their employment in specific problems is far from being trivial. In this Letter I describe a new formulation for 3D optical field quantization, and discuss its application to optical parametric amplification.

All the machinery is based on the paraxial and slowly varying envelope approximation (SVEA), with reference to vacuum and to a normally dispersive medium. However a general theory can be envisaged. The basic element is recognizing that a traveling pulsed beam can be expressed as a continuous set of harmonic oscillators, each weighted by a rigidly moving 3D wave-packet. It corresponds to a non-dispersive non-diffracting solution of the relevant classical wave equation, with velocity $`v`$. This is the so-called X-wave transform, or X-wave expansion, as outlined below. Lu00 ; Salo01b  Then the quantization comes naturally, following standard approach, using continuous-mode operators. DiracBook ; Blow90

X-waves are 3D invariant pulsed-beams traveling without diffraction and dispersion, also known to as “progressive undistorted waves” (see e.g. the recent reviews Saari01 ; Recami03  and references therein). Their appearance in nonlinear optical processes has been experimentally investigated in DiTrapani03 .

For free propagation, it turns out that the quasi-monochromatic beam, with carrier angular frequency $`\omega `$, behaves like a quantum gas of particles which, in vacuo, have mass $`m`$ satisfying the relation $`mc^2=\mathrm{}\omega `$. The rigidly moving “eigenmodes” (the X-waves), while not having a finite norm (like plane waves), do have a strong degree of localization, thus addressing the issue of spatio-temporal photon confinement (see e.g. Chan02  and references therein). In a simple way, this approach provides the quantum description of the 3D evolution of a classical, paraxial and slowly varying, beam. Furthermore classical infinite-energy X-waves, here described as coherent states, do have finite expectation values for energy in the quantized picture.

## II Derivation of the paraxial equations

In this article a new quantum-mechanical description of pulsed light beams is introduced. In order to avoid as much as possible unnecessary complexities and involved formalism, paraxial slowly varying beams are considered. Given the fact that progressive undistorted waves are well-studied solutions of the wave equation, a much more general theory can be envisaged. Typical experiments in nonlinear optics are described in the framework of the considered approximations and, for similar reasons, the treatment is limited to circularly symmetric beams.

The basic model can be obtained by the scalar 3D wave equation for the electromagnetic field $`E`$ in a medium with refractive index $`n`$:

$$^2E\frac{n^2}{c^2}\frac{^2E}{t^2}=0\text{.}$$
(1)
By letting $`E=Re[A\mathrm{exp}(i\omega t+ikz)]/\sqrt{ϵ_0n^2/2}`$, with $`k=\omega n/c`$, in the paraxial approximation the usual Foch-Leontovich equation can be derived:

$$2ik\frac{A}{z}+2i\frac{\omega }{c^2}\frac{A}{t}+_{}^2A\frac{1}{c^2}\frac{^2A}{t^2}=0\text{.}$$
(2)
Equations like (2) are typically adopted in nonlinear optics, and interpreted as a propagation problem along $`z`$. If $`A`$, as a function of $`t`$, is slowly varying, they can be equivalently cast as a time-evolution problem; the latter is the formulation adopted here in order to be as much as possible consistent with standard quantum mechanics. With this aim, introducing $`\zeta =zct`$ and $`\tau =t`$, Eq. (1) is rewritten as

$$2i\frac{\omega }{c^2}\frac{A}{\tau }+_{}^2A\frac{^2A}{\zeta ^2}=\frac{1}{c^2}\frac{^2A}{\tau ^2}\frac{2}{c}\frac{^2A}{\tau \zeta }\text{.}$$
(3)
If the SVEA is applicable, the right-hand side of the previous equation is negligible and, coming back to the original variables $`(z,t)`$, Eq. (3) can be written as

$$i\frac{A}{t}+i\omega ^{}\frac{A}{z}\frac{\omega ^{\prime \prime }}{2}\frac{A}{z^2}+\frac{\omega ^{}}{2k}_{}^2A=0\text{,}$$
(4)
while being $`\omega ^{}=c`$ and $`\omega ^{\prime \prime }=c^2/\omega `$. Eq (4) can be treated as an evolution problem with respect to time $`t`$.

The very same procedure can be repeated in a dispersive medium, as done in Ref. Kennedy88, , while retaining terms up to the second order dispersion. The model (4) still holds with $`\omega ^{}`$ and $`\omega ^{\prime \prime }`$ the derivatives of the dispersion relation $`\omega =\omega (k)`$. Hence, since it can describe vacuum as well as a dispersive material, 111For uniformity, I will limit the treatment to a normally dispersive material ($`\omega ^{\prime \prime }>0`$). and it is typically used with reference to experiments with pulsed-diffracting beams, it will be the starting point for the quantization procedure. Without loss of generality the envelope field $`A`$ has been normalized so that the time averaged electromagnetic energy in the whole 3D space is

$$=|A|^2𝑑x𝑑y𝑑z\text{,}$$
(5)
i.e. any type of box-normalization is not necessary.

## III The X-wave transform and X-wave expansion

Assuming radial symmetry hereafter ($`r\sqrt{x^2+y^2}`$), the general solution for the initial value problem is expressed by the Fourier-Bessel integral

$$A=_{\mathrm{}}^{\mathrm{}}_0^{\mathrm{}}k_{}J_0(k_{}r)S(k_{},k_z)e^{ik_zZi\mathrm{\Omega }t}𝑑k_{}𝑑k_z\text{,}$$
(6)
with $`Z=z\omega ^{}t`$. Eq. (6) furnishes the field at instant $`t`$, given its spectrum $`S`$ at $`t=0`$, with transversal and longitudinal wave-numbers $`k_{}`$ and $`k_z`$, and $`\mathrm{\Omega }=\omega ^{\prime \prime }k_z^2/2+\omega ^{}k_{}^2/2k`$. With the following change of variables ($`\alpha >0`$)

$$\begin{array}{c}k_{}=\alpha \sqrt{\omega ^{\prime \prime }k/\omega ^{}}\hfill \\ k_z=\alpha v/\omega ^{\prime \prime }\text{,}\hfill \end{array}$$
(7)
I have, after some manipulations,

$$A=_{\mathrm{}}^{\mathrm{}}e^{i\frac{v^2}{2\omega ^{\prime \prime }}t}\mathrm{\Psi }_v(r,Zvt)𝑑v\text{,}$$
(8)
while being

$$\mathrm{\Psi }_v(r,Zvt)=_0^{\mathrm{}}X(\alpha ,v)J_0(\sqrt{\frac{\omega ^{\prime \prime }k}{\omega ^{}}}\alpha r)e^{i(\alpha \frac{v}{\omega ^{\prime \prime }})(Zvt)}𝑑\alpha \text{.}$$
(9)
Hence $`A`$ can be expressed as a superposition of self-invariant beams (i.e. waves traveling along the $`z`$direction without dispersion and diffraction) with different velocities. The quantity $`X(\alpha ,v)=\frac{k\alpha }{\omega ^{}}S(\alpha \sqrt{\frac{\omega ^{\prime \prime }k}{\omega ^{}}},\alpha \frac{v}{\omega ^{\prime \prime }})`$ corresponds to the so-called X-wave transform.Lu00  X-waves are defined as the propagation invariant solutions of the wave equation (4), for more details see Ref. Conti03,  and references therein.

$`X(\alpha ,v)`$ can be expanded by Laguerre functions Salo01b  $`L_p^{(1)}(2x)e^x`$ with respect to $`\alpha `$:

$$\begin{array}{c}X(\alpha ,v)=\mathrm{\Sigma }_{p=0}^{\mathrm{}}C_p(v)f_p(\alpha )\hfill \\ \\ f_p(\alpha )=\sqrt{\frac{k}{\pi ^2\omega ^{}(p+1)}}(\alpha \mathrm{\Delta })L_p^{(1)}(2\alpha \mathrm{\Delta })e^{\alpha \mathrm{\Delta }}\text{,}\hfill \end{array}$$
(10)
where $`\mathrm{\Delta }`$ is a reference length that is related to the spatial extension of the beam (the specific form of $`f_p`$ is chosen for later convenience). $`A`$ is written as a superposition of basis orthogonal X-waves (given by (9) with $`X(\alpha ,v)=f_p(\alpha )`$):

$$A=\mathrm{\Sigma }_{p=0}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}C_p(v)e^{i\frac{v^2}{2\omega ^{\prime \prime }}t}\psi _p^v(r,Zvt)𝑑v\text{.}$$
(11)
The following relation holds by construction

$$<\psi _q^u(r,Zut)|\psi _p^v(r,Zvt)>=\delta _{pq}\delta (uv)\text{,}$$
(12)
where $`<f|g>`$ denotes the integral over $`x,y,z`$ of $`f^{}g`$. Eq. (12) reveals that the X-waves have infinite norm, like the plane waves typically adopted for field quantization.

## IV Quantization

The classical energy of the pulsed beam, Eq. (5), is given by

$$=\mathrm{\Sigma }_{p=0}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}|C_p(v)|^2𝑑v\text{,}$$
(13)
while the field in (11) can be rewritten as

$$A=\mathrm{\Sigma }_{p=0}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}C_p(v,t)\psi _p^v(r,Zvt)𝑑v\text{,}$$
(14)
being

$$\frac{dC_p}{dt}(v,t)=i\omega _p(v)C_p(v,t)$$
(15)
with $`\omega _p=v^2/2\omega ^{\prime \prime }`$. Hence $`A`$ is the integral sum of harmonic oscillators, with complex amplitudes $`C_p`$, each associated to a traveling invariant beam (corresponding to the usual cavity mode, with the difference that it is rigidly moving). Introducing the real-valued “position” $`Q_p(v)`$ and “momentum” $`P_p(v)`$, by letting

$$C_p(v,t)=\frac{1}{\sqrt{2}}[\omega _p(v)Q_p(v,t)+iP_p(v,t)]\text{,}$$
(16)
$``$ is written as the sum of the time-independent classical energies of each oscillator

$$=\mathrm{\Sigma }_{p=0}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\frac{1}{2}[P_p(v,t)^2+\omega _p(v)^2Q_p(v,t)^2]𝑑v\text{,}$$
(17)
$`\omega _p(v)`$ is the resonant frequency of each “mode”.

Following standard approach to quantization, ItzyksonBook ; LoudonBook  it is natural to write the Hamiltonian like (omitting extrema for the integral and the $`\mathrm{\Sigma }`$ symbol)

$$H=\mathrm{\Sigma }_p\mathrm{}\omega _p(v)a_p^{}(v)a_p(v)𝑑v=\mathrm{\Sigma }_p\frac{mv^2}{2}a_p^{}(v)a_p(v)𝑑v\text{,}$$
(18)
with $`[a_q(u),a_p^{}(v)]=\delta _{pq}\delta (uv)`$. Note that the zero-point energy has been omitted in (18), as a standard renormalization in quantum field theory. ItzyksonBook

Some comments are needed at this point. First of all I observe that (18) corresponds to a set of freely moving particles, each contributing to the energy by the kinetic term, with mass defined by

$$m=\frac{\mathrm{}}{\omega ^{\prime \prime }}\text{.}$$
(19)
Thus the freely propagating pulsed beam is represented like a quantum fluid, or gas.Landau  The eigenstates are given by $`|l,n,v>=[a_l^{}(v)]^n|0>`$, denoting $`n`$ particles with velocity $`v`$ in the traveling “mode” (or fundamental X-wave) $`\psi _l^v`$, and energy $`nmv^2/2`$. Then, it is to note that, differently from standard heuristic approaches to field quantization, the Hamiltonian directly corresponds to the energy of the beam, and not to the second quantized form, i.e. to the Hamiltonian functional giving Eq. (4).

The optical field amplitude operator is given by

$$A=\mathrm{\Sigma }_pe^{\frac{i}{\mathrm{}}\frac{mv^2}{2}t}\psi _p^v[z(\omega ^{}+v)t,r]\sqrt{\mathrm{}\omega _p(v)}a_p(v)𝑑v\text{.}$$
(20)
The mean value of $`A`$ and of the energy density, in the state with one elementary excitation on the $`l`$fundamental X-wave at velocity $`v`$, are: $`<l,1,v|A|l,1,v>=0`$ and $`<l,1,v|A^{}A|l,1,v>=\mathrm{}\omega _l(v)|\psi _l^v(Zvt,r)|^2`$.

Coherent states can we constructed as well. For example, taking (with obvious notation)

$$|\alpha _l(v)>=e^{\frac{|\alpha _l(v)|^2}{2}}\mathrm{\Sigma }_{n=0}^{\mathrm{}}\frac{\alpha _l(v)^n}{\sqrt{n!}}|l,n,v>\text{,}$$
(21)
with $`\alpha _l(v)`$ an arbitrary complex number, I have a state that has an expectation value for the optical field profile given by the fundamental X-wave $`\psi _l^v(Zvt,r)`$ and corresponds to a classical X-wave. 222A generic optical field will be given by a superposition of these states. However there is a remarkable difference: the expectation value of the energy is $`\mathrm{}\omega _l(v)|\alpha _l(v)|^2`$. While this state has infinite energy in a classical context, in a quantum mechanical picture X-waves do have finite energy.

On the other hand it is well known, from quantum field theory, that $`|l,n,v>`$ is not normalizable, and $`a_l^{}(v)`$ must be intended as operator valued distributions. ItzyksonBook  This is due to the fact that a state with a definite $`v`$ is an idealization, like an elementary particle with a definite momentum $`k`$. Any experimental apparatus implies a spread, and an indeterminacy in the velocity (corresponding to a normalized state obtained as a continuous superposition of $`|l,n,v>`$).

By particularizing Eq. (19) to the case of vacuum, using the previously mentioned relations, it is found

$$mc^2=\mathrm{}\omega \text{,}$$
(22)
which shows that the effective mass of the quantum gas particles is determined by the carrier frequency of the beam via the Einstein relation.

It is also notable that, in vacuo ($`\omega ^{}=c`$), superluminal components are encompassed in Eq. (20). But this is not surprising. Indeed the superluminal propagation of X-waves does not lead to paradoxes, since they are stationary solutions (i.e. they fill all the space, like plane waves) of the wave equations (see also Ciattoni03 ; Recami03 ). This can be shown to be perfectly consistent with the principles of standard special relativity, in the framework of the covariant Maxwell equations.CiattoniPrivate  Furthermore it is well known that superluminal propagation has a not zero probability, and hence it is predicted, in standard quantum electrodynamics. Feynman  Single-photon superluminal effects have been experimentally investigated in Ref. Steinberg93, . 333Obviously the paraxial-SVEA model is not Lorentz invariant, the relevant generalization will be the subject of forthcoming papers.  Before considering an application of this formalism I want to underline that a completely different wavelet quantization of relativistic fields has been previously developed by Kaiser.Kaiser03

## V Optical Parametric Amplification

I consider the phase-matched co-directional optical parametric amplification (OPA) Hong85  of two frequencies, $`\omega _1`$ and $`\omega _2`$, attainable in quadratic, or cubic, nonlinear media. For simplicity the pump beam is treated as a constant term in the Hamiltonian. A quantized field is associated to $`\omega _1`$ and $`\omega _2`$, with parameters $`\omega _{1,2}^{}`$, $`k_{1,2}=\omega _{1,2}n_{1,2}/c`$. $`\omega ^{\prime \prime }`$ is taken equal for both of them, and $`\omega _1^{}\omega _2^{}`$. $`A_{1,2}`$ denote the corresponding optical quantum fields:

$$\begin{array}{c}A_1=\mathrm{\Sigma }_pe^{i\frac{v^2t}{2\omega ^{\prime \prime }}}\psi _{p,1}^v[z(\omega _1^{}+v)t,r]\sqrt{\mathrm{}\omega _p(v)}a_p(v)𝑑v\hfill \\ A_2=\mathrm{\Sigma }_pe^{i\frac{v^2t}{2\omega ^{\prime \prime }}}\psi _{p,2}^v[z(\omega _2^{}+v)t,r]\sqrt{\mathrm{}\omega _p(v)}b_p(v)𝑑v\text{.}\hfill \end{array}$$
(23)
The basis X-waves are denoted by $`\psi _{p,j}^v`$ in order to indicate that they have to be calculated after Eq. (10) with parameters ($`\omega _j^{}`$ and $`k_j`$) of $`\omega _j`$; the corresponding spectra are denoted by $`f_{p,j}(\alpha )`$ ($`j=1,2`$).

The total Hamiltonian is given by the sum of the free terms and the interaction Hamiltonian $`H_I`$:

$$H=\mathrm{\Sigma }_p\mathrm{}\omega _p(v)(a_p^{}a_p+b_p^{}b_p)𝑑v+H_I\text{.}$$
(24)
$`H_I`$ can be determined by the classical counterpart given by $`_I=\chi <A_1|A_2^{}>+c.c.`$ ($`c.c.`$ denotes complex conjugation). After some (lengthy but straightforward) manipulations ($`h.c.`$ is hermitian conjugate):

$$H_I=\mathrm{}\mathrm{\Sigma }_{pq}𝑑u𝑑v\chi _{pq}(u+v)\sqrt{\omega _p(v)\omega _q(u)}e^{iF(u,v)t}a_p^{}(u)b_q^{}(v)+\text{h.c.}$$
(25)
where $`F(u,v)=(u^2+v^2)/2\omega ^{\prime \prime }+(uv+\omega _1^{}\omega _2^{})(v\rho u)/\omega ^{\prime \prime }`$, with $`\rho =(k_1\omega _2^{}/k_2\omega _1^{})^{1/2}`$. The interaction term $`\chi _{pq}(\nu )`$ is

$$\chi _{pq}(\nu )=\frac{4\pi ^2\chi }{\nu }\sqrt{\frac{\omega _1^{}\omega _2^{}}{k_1k_2}}f_{p,1}[\frac{\nu }{(1+\rho )\omega ^{\prime \prime }}]f_{q,2}[\frac{\rho \nu }{(1+\rho )\omega ^{\prime \prime }}]\theta (\nu )\text{,}$$
(26)
with $`\theta (\nu )`$ the unit step function.

By applying standard time-dependent perturbation theory, assuming as the initial state $`|0>`$ at $`t=0`$, at first order it is

$$\begin{array}{c}|\psi ^{(1)}(t)>=\mathrm{\Sigma }_{pq}e^{iK(u,v)t}\chi _{pq}(u+v)G(u,v,t)\sqrt{\omega _p(v)\omega _q(u)}a_p^{}(u)b_q^{}(v)𝑑u𝑑v|0>\text{,}\hfill \end{array}$$
(27)
with

$$\begin{array}{c}K(u,v)=\frac{(uv+\omega _1^{}\omega _2^{})(v\rho u)}{2(1+\rho )\omega ^{\prime \prime }}\hfill \\ g(u,v)=\frac{(u^2+v^2)}{\omega ^{\prime \prime }}+\frac{(uv+\omega _1^{}\omega _2^{})(v\rho u)}{(1+\rho )\omega ^{\prime \prime }}\hfill \\ G(u,v,t)=\frac{2\mathrm{sin}[g(u,v)t/2]}{g(u,v)}\text{.}\hfill \end{array}$$
(28)
The state given in Eq. (27) is a continuous-variables entangled superposition of particles at the two generated frequencies, traveling with different velocities. Indeed the weight function is not separable with respect to $`u`$ and $`v`$ (see e.g. Ref. Chan02,  and references therein).

The transition probability to the two particles state $`a_p^{}(u)b_q^{}(v)|0>`$ is given by

$$𝒫_{pq}(t,u,v)=\omega _p(v)\omega _q(u)|\chi _{pq}(u+v)|^2\frac{sin[g(u,v)t/2]^2}{[g(u,v)t/2]^2}t^2.$$
(29)
As $`t\mathrm{}`$, $`𝒫_{pq}`$ tends to a Dirac-delta peaked at $`g(u,v)=0`$. Entangled particles, with velocities $`u`$ and $`v`$, are associated, for large propagation distances, to a point in the plane $`(u,v)`$, lying on the parabola $`g(u,v)=0`$. 444Related results, at a classical level, seem to be those of Ref. Picozzi02, , where coherence along space-time trajectories in parametric wave mixing has been investigated. Of particular relevance is the low velocity region.

Proceeding as in quantum fluid theory, the small momenta approximation can be applied.Landau  This means neglecting the quadratic terms in $`u,v`$ in $`g(u,v)`$, which corresponds to small velocities, consistently with the SVEA. In this case $`g(v\rho u)(\omega _1^{}\omega _2^{})/(1+\rho )\omega ^{\prime \prime }`$, showing that, as $`t\mathrm{}`$,

$$𝒫_{pq}(t,u,v)\omega _p(v)\omega _q(u)|\chi _{pq}(u+v)|^2\frac{2\pi (1+\rho )\omega ^{\prime \prime }t}{|\omega _1^{}\omega _2^{}|}\delta (v\rho u)\text{.}$$
(30)
Thus, in the asymptotic state the two particles travel approximately at the same velocity ($`\rho 1`$ in practical cases).

For large propagation distances, the quantum fluid is hence a superposition of entangled pairs locked together. This gives a clear picture of what entanglement is to be intended for. After the generation, the two excitations are associated to X-shaped spatial distributions of energy that travel at the same $`v`$. Since the velocity determines the angular aperture (far from the origin) of the double cone which forms the X-wave,Ciattoni03  a three-dimensional region where the two particles can be revealed simultaneously is determined. The duration of the temporal profile fixes the useful time window. In some sense, the paired elementary excitations are delocalized in an highly localized 3D X-shaped rigidly moving region.

## VI Conclusions

Using X-waves as basis, a traveling wave quantization leads to the representation of a pulsed beam as a quantum gas of free particles. A structured, localized, spatio-temporal distribution of energy is associated to each fundamental excitation. These elements furnish an alternative point of view regarding the concept of “photon”, and its localization, without introducing box-normalization or related ideas.

A 3D quantum propagation is reduced to a one dimensional evolution, thus enabling simple treatments of nonlinear optical processes. This is mainly a consequence of the paraxial approximation, which unavoidably fixes a preferential direction. For this reason, this approach is well suited for describing optical experiments. It is also susceptible of generalizations in Lorentz covariant theories.

During optical parametric amplification the particles composing the pulsed beam become entangled in pairs moving, approximately, at the same velocity. These excitations are delocalized into a 3D rigidly moving conical region, where they can be revealed simultaneously. The role of first order and of second order dispersion is clearly stated, while including diffraction. This provide a complete spatio-temporal picture of a nonlinear quantum process.

## Aknowledgements

I thank G. Assanto, A. Ciattoni, B. Crosignani, and S. Trillo for invaluable discussions.


format: text

text:
# Autonomous Energy Transducer: Proposition, Example, Basic Characteristics

## I Introduction

How energy is converted from one form to another form has been extensively studied in these days. In biochemical reaction, some molecules are known to work for “energy transducer”. Such energy transduction is important in enzyme function, molecular motor, and so forth, while recent progress in nanotechnology makes us possible even to design a microscopic machine that works as such energy transducer.

The question whether such molecular machine functions in the same way as a macroscopic machine was first addressed up by Oosawa  as a dichotomy between tight and loose coupling between input and output. In the tight coupling, chemical energy is converted into mechanical work via definitely scheduled successive conformational changes of molecules, and at each stage the conformational change has one-to-one correspondence to a mechanical work. This mechanism is similar to most of our macroscopic machine. Consider a typical machine with gears. It converts some form of energy to work via several rigid gears, following external operations. Each gear transforms one kind of motion to another directly and instantaneously, where fluctuations in the transformation process are negligible.

In contrast, in the loose coupling mechanism, the manner or the route to establish energy conversion does fluctuate by each event, and the energy is not necessarily converted all in once at each stage. Here, the amount of output work is not precisely specified, but varies by each event. By repeating trials, the output work extracted from the same amount of input energy (e.g., by a single reaction from ATP to ADP) is distributed with some fluctuations. Oosawa proposed that this loose coupling is a general property of a system that robustly works under large thermal fluctuations of the configuration of molecules even though the input energy is of the same order with the energy of thermal fluctuations.

This problem is first posed in the study of molecular motors, where conversion from chemical energy to directional mechanical work within a fluctuating environment is studied. Indeed, there are some experimental results suggesting loose coupling mechanisms. For instance, the energy is stored in a molecule over an anomalously long time between the chemical reaction event and emergence of directional motion , and multiple steps of conversion to work are repeated during molecular processes within one cycle of the chemical reaction . There the resulting output steps is widely distributed as is expected in loose coupling mechanism.

Whether the energy conversion of a biomotor is loose or tight may depend on detailed molecular structure. Still, the question how the energy conversion in loose coupling mechanism works is quite general, and can be discussed in any interacting molecules. One possible mechanism for the loose coupling was proposed as thermal ratchet at a statistical level , although dynamical process for the loose energy conversion is not yet clarified. Little is known what types of loose conversion are possible, or in what condition and by what mechanism the loose coupling functions.

In the present paper, we propose a scheme of dynamical energy transduction of autonomous regulation in highly fluctuating configurations, by introducing a system consisting of interacting molecules with internal dynamics. The specific choice of a model we adopt for this energy conversion is inspired by some experiments suggesting oscillatory molecular dynamics at a time scale of molecular work , whereas we do not pursue to make a direct correspondence with the experiment. Rather, we adopt the model, expecting that the proposed mechanism itself is rather general. By setting up an example of a dynamical systems model that works as a loose coupling, we show that energy transduction is possible within fluctuating environment. By analyzing the dynamical process of this energy transduction, we find some characteristic features of the dynamics that makes the coupling in this scheme generally loose. Robust transduction of energy within fluctuating configuration is discussed from the viewpoint of dynamical systems. The proposed mechanism may also shed new light on the function of enzyme in chemical reaction, design of a nano-machine that works at a molecular scale, and biomotors.

## II Definition of Autonomous Energy Transducer at Molecular Scale

To study this problem, we start by defining a machine that works at a molecular scale, with which we are concerned here. The requisite here is essential to a system that works with large fluctuations of configurations in molecules.

Consider a system composed of several degrees of freedom. ‘Input’ is injected into the system through some restricted degrees of freedom of it, that are pre-defined. Then, as a reaction to the input, there arises dynamical change in the system, which leads to ‘output’ as a typical motion of another set of restricted degrees of freedom. Here, the degrees of freedom for the input and for the output are distinct. Statistically, the output brings directional flow or mechanical work, although individual events of the output would be rather probabilistic. The dynamics of the system is given by the interaction among the degrees of the freedom, as well as noise from the heat bath. Through the dynamics, input energy is transformed to (stochastic) output motion, that gives an output work statistically. This energy transducer works with large fluctuations of configurations of the system, without any external control, once the energy is injected as input. Accordingly we postulate the following properties:

* Details in input event are not controlled. Neither the timing of the injection with regards to the molecule fluctuations, nor the direction of the input motion is specified precisely. Even the magnitude of the input is not restricted.
* Once the input is injected, no external control is allowed, and the process is governed by a given rule for temporal evolution, while fluctuations at a molecular scale are inevitable.
* In spite of this uncontrolled condition, mechanical work or directional flow can be extracted from the output motion, on the average. The system adapts itself to produce output, against the change of the input timing, direction, or magnitude.

Then, we address the question how such system is possible. In the present paper, we provide a specific example of a system with such loose-coupling, and discuss its characteristic features which realize the above three conditions. The characteristic features of dynamical process are summarized as follows:

After the energy is given as the input, the input energy is sustained with very weak dissipation for some time, and is not directly transduced to the output. This suppression of energy flow into other modes comes from some restriction to the interaction among the system. On the one hand, this restriction makes a change of molecular dynamics from the fluctuating motion to the directional motion later.

During the time while the energy is stored, detailed information of the input (such as the direction or the timing), except the energy value, is lost due to noisy chaotic dynamics. This makes the output distributed. In other words, the input-output relationship is loose, in the sense that the output (i.e., displacement or extracted work) is distributed around its mean. There is a rather general relationship between the variance and the mean of the outputs, as will be shown.

By studying a dynamical model with several degrees of freedom, we have confirmed the above features, when the model satisfies the requisite (i)-(iii), to extract mechanical work from the energy flow within highly fluctuating configurations. On the other hand, no substantial work is extracted from the output, when these features are not observed in the models. These results imply that the autonomous energy transduction works as a loose coupling mechanism.

## III Model

As a specific system to satisfy the condition (i)-(iii), we consider the following model ***The present model is originally inspired by the study of motor protein (e.g., myosin) that moves along the chain of the rail protein, actin. . The motor part consists of a ‘head’, of position $`x_h`$, and an internal ‘pendulum’ represented by $`\theta `$. The motor interacts with rail, represented by a one-dimensional chain of lattice sites positioned at $`x_i`$ with index $`i`$. The internal pendulum is excited when energy is applied to the system. Every degree of freedom, except for $`\theta `$, is in contact with a heat bath, generating random fluctuations described by a Langevin equation with damping. The interaction potential $`V(x_hx_i,\theta )`$ between the chain and the head is spatially asymmetric and its form depends on the angle of the pendulum (see Fig.1). The periodic lattice is adopted as the chain, to study directional motion in an asymmetric periodic potential, as is often studied in the study of thermal ratchets. Every degree of freedom, except for the internal pendulum, is in contact with the heat bath, generating random fluctuations described by a Langevin equation.

The equations of motion for this system are chosen as

$`m_c\ddot{x}_i`$ $`=`$ $`\gamma \dot{x}_i+\sqrt{2\gamma T}\xi _i(t)`$ (1)
$``$ $`K_c\left\{(x_iiL)+(x_i{\displaystyle \frac{x_{i1}+x_{i+1}}{2}})\right\}{\displaystyle \frac{V(x_hx_i,\theta )}{x_i}},`$ (2)
$`m_h\ddot{x}_h`$ $`=`$ $`\gamma \dot{x}_h+\sqrt{2\gamma T}\xi _h(t){\displaystyle \underset{i}{}}{\displaystyle \frac{V(x_hx_i,\theta )}{x_h}},`$ (3)
$`m_\theta \ddot{\theta }`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{V(x_hx_i,\theta )}{\theta }},`$ (4)
where $`T`$ is the temperature, $`\gamma `$ is a friction coefficient and $`\xi _\alpha (t)`$ represents Gaussian white noise. Here, we use the units Boltzmann constant $`k_B=1`$. $`K_c`$ and $`L`$ are the spring constant and the natural interval between two neighboring lattice sites in the chain. To observe directional motion of the head, the chain is also connected to a fixed ground via a spring with a constant $`K_c`$. Here, $`m_c`$, $`m_h`$ and $`m_\theta `$ are mass of the respective degrees of freedom, $`m_c=m_h=1`$ and $`m_\theta =O(0.01)`$. The inertia is not ignored, as small friction coefficient $`\gamma `$ is adopted. Although this could be reasonable for the application to enzymatic reaction or to a design of nano-machine, one might think that biomotors should be treated as overdamped systems, considering that for small object that moves with a slow time scale, water is highly viscous fluid. Here, we introduce inertia term, because several experimental results by Yanagida’s group may not necessarily match with this standard view . Also, each variable in our model does not necessarily correspond to the atomic configuration of a protein. Rather it is a mesoscopic variable, for which effective inertia may be relevant to long-term storage of energy in dynamics of proteins, and allow for oscillatory dynamics among interacting degrees of freedom. We will come back to this ’inertia problem’ in §8 again.

The potential form is asymmetric in space as shown in Fig.1, where the characteristic decay length of the interaction is set at a smaller value than $`L`$, so that the interaction is confined mostly to the nearest lattice sites. Here we adopt the following potential form, $`V(\mathrm{\Delta }x,\theta )=V_x+V_\theta `$ with $`V_x=K_h\mathrm{tanh}(p(\mathrm{\Delta }xr))/\mathrm{cosh}(d\mathrm{\Delta }x),`$ and $`V_\theta =K_h(1\mathrm{cos}\theta )/2\mathrm{cosh}(d\mathrm{\Delta }x).`$ Here the parameters $`p`$ and $`r`$ determine the degree of asymmetry, while $`K_h`$ and $`d^1`$ give the coupling strength and decay length of the interaction, respectively. Specific choice of this form is not important. We have simulated our model choosing several other potential forms with asymmetry, for instance $`V(\mathrm{\Delta }x,\theta )=K_h(1\mathrm{cos}\theta )/(\mathrm{exp}(p\mathrm{\Delta }x)+\mathrm{exp}(d\mathrm{\Delta }x)),`$ and obtained qualitatively the same results with regards to the directional motion. In this paper, the parameters for the potential are fixed at $`p=10`$, $`r=0.3`$, $`K_h=0.2`$ and $`d^1=0.25`$.

The pendulum mode without direct coupling to the heat bath is adopted here just as one representation of long-term storage of energy  experimentally suggested in molecular motors . As long as there exists such slow relaxation mode, our results follow, and this specific choice is just one example for it. Any mode realizing slow relaxation can be adopted instead.

## IV Output directional motion

Once a sufficient energy is injected to the pendulum to drive the system out of equilibrium, we have observed that the head moves a few steps in one direction, on its way in the relaxation to thermal equilibrium . An example of the time course of the head motion is given in Fig.2. How many steps the head moves after the energy injection depends on each run with different seed for the random force, even though the values of parameters including $`E_0`$ are kept constant.

In Fig. 3(a), the average step size is plotted as a function of temperature ($`1/T`$). The mean step size only weakly depends on the temperature for small $`T`$. It shows only slight decrease by decreasing the temperature. This is in strong contrast with the thermally activated process, in which the rate of crossing the potential barrier (for the head to shift) is given by $`\mathrm{exp}(U_0/T)`$ with an activation energy $`U_0`$, according to Arrhenius law. At thermal equilibrium, indeed, the head shows such thermally activated motion. The resulting diffusion rate of the head Brownian motion obeys this exponential form as shown in Fig.3(b), whereas the values of $`U_0`$ agree with the activation energy derived analytically .

The dependence in Fig.3(a) is distinct from such dependence, and the motion of the head is not thermally activated. This is natural, since the energy for the directional motion is extracted from the injected energy $`E_0`$. As a result, the head motion possesses different properties from the thermally activated motion as follows. For a fixed temperature, the average step of the head increases monotonically (almost linearly) with $`E_0`$, when the injected energy $`E_0`$ is larger than a threshold (see Fig.4(a)). This threshold is not equal to $`U_0`$ for the Brownian motion, but shows the following dependence upon temperature: For higher temperature, this threshold energy goes to zero, and any small injected energy can lead to the directional motion on the average, whereas, with decreasing temperature, the threshold energy tends to converge to a certain value around $`0.3`$. This result suggests that there is a certain minimum energy required for the head to shift to the next lattice site.

The temperature dependence suggests that the head motion is not thermally activated. Rather, it seems that the motion takes place when a suitable temporal correlation exists between the head motion and thermal fluctuations of the chain. It should also be noted that the directional motion is clearly observed even by the injection of the energy of $`O(10)`$ times of $`k_BT`$. For example, the average motion is about 2 steps, for the injection of the energy of $`20k_BT`$, for $`T=0.02`$.

The linear increase of the step size with $`E_0`$ implies that the injected energy is not wasted even if it is larger than the required threshold value for the directional motion of the head. This implies that the injected energy is converted to the directional motion not at once but step by step with the successive step motion of the head. This feature is revisited in Section 6.

To make this successive energy conversion possible, it is important that the chain be sufficiently flexible. Indeed, as is shown in Fig.5, if the spring constant $`K_c`$ is too large, directional motion of the head is not observed. The directional motion is not possible, if the chain is too tight, or too flexible. The average step size has a peak around $`K_c0.5`$. This suggests that a flexible motion that appears under a proper stiffness, plays an important role to produce the directional motion as an output of the energy transduction. Indeed, around this parameter value, there is a transition in the dynamics of the head and neighboring few lattice points, to show correlated motion, and accordingly a rather long-range correlation appears there.

## V Characteristics of Dynamic Energy Transduction with Loose Coupling

Now we discuss how the above directional motion is achieved in the term of dynamical systems. We note that the present mechanism works only as a loose coupling mechanism in this section.

First, let us discuss how the system ‘forgets’ initial conditions and how broad distribution of output step motion appears. Before the directional motion of the head, the motion of the triplet {head, pendulum, rail } shows a weakly chaotic motion, with which the input information (e.g., the configuration of the system when the input is applied) is lost. The chaotic motion here only works as weak perturbation to the two-body problems consisting of the interactions between the head and the pendulum and that between the head and the adjacent lattice site. On the other hand, when the head exhibits a step motion, there is strong interaction between the head and the lattice sites, that lead to a stronger chaotic instability arising from three-body or five-body (including the neighboring lattice points) motions .

The change of the degree of instability is demonstrated in Fig.6, where the evolution of difference between two close orbits is plotted. As shown, there are two rates of exponential divergence, one that continues during the energy storage, and the other with the flow of energy from the pendulum part. The former is given by perturbations to the two-body motions, and the latter is induced by strong interaction with the three- or five-body motion. The orbital instability here is examined by taking two identical initial configurations, and applying identical input energy and identical random noise (i.e., sequence of random number). A slight disturbance is added to the reference system at the timing of energy injection. The evolution of these two systems with tiny difference in configurations results in quite different number of output steps (see Fig.6(a)). Then, it is impossible to predict or control the number and the direction of steps by the configuration of the system at the excitation. On the other hand, the directional motion appears rather independently of the configuration and, therefore the energy conversion of the system is robust against thermal noise, and disturbances in the input event.

One consequence of these chaotic dynamics together with the influence of the heat bath is a loose relationship between the input and the output. In fact, the output number of steps shows a rather broad distribution, as shown in Fig.7, where one step means the displacement of the head for one lattice site along the chain. The distribution was obtained by taking $`1000`$ samples with arbitrarily chosen configurations at the moment of energy injection, while satisfying thermal equilibrium, and also by taking different random sequences $`\xi _\alpha (t)`$. The resulting distribution is rather broad, although the input energy $`E_0`$ is identical for all the samples. The difference in the output step comes mainly from intrinsic chaotic dynamics of the system mentioned above.

In the distribution of the number of steps, we note that at low temperature the mean displacement is almost equal to its variance (see Fig.7(b)). If the step motion is a result of rare stochastic events, one could expect Poisson distribution. Indeed, if one disregards the small tail of the distribution at negative steps, the obtained distribution is not far from Poissonian, for low temperature. Even for high temperature, the linear relationship between the mean and variance is still valid with an offset depending on the temperature. This Poisson-type distribution suggests that random process underlies the output head motion, as expected from the above chaotic dynamics.

The relationship in Fig.7(b) indicates that to obtain sufficient output amount in directional motion, large fluctuation is inevitable. Indeed, it is not possible, at least in the present model, to make the step distribution sharp while keeping sufficient directional motion of the head, even if the parameter values are changed. For example, one could make our model rigid, by increasing stiffness $`K_c`$ in the rail and suppressing the fluctuations of the configuration. Then, the directional motion is highly suppressed, as already shown in Fig.5. The average displacement as a function of the stiffness has a maximum around $`K_c0.5`$, where the fluctuation of the displacement of rail elements is rather large. On the other hand, for $`K_c=1.3`$, i.e., for a rigid spring of the rail, the average displacement is only 2% of the maximum. There, most of the input energy then is wasted. Accordingly, to achieve larger directional steps without wasting input energy, the step distribution has to be broad. This demonstrates the relevance of loose coupling mechanism to energy transduction.

This relevance of loose coupling is expected to be general under our postulates for the molecular transducer. Note that rigid process in tight coupling mechanism can transfer energy very well when configurations of every two parts between which the energy is transferred is precisely determined, but the efficiency will go down drastically if this precise condition for the configurations is not satisfied. In the molecular energy transducer we are concerned, such precise control of configurations of molecules is impossible. Furthermore, as the detailed information on the input is lost, control of the motion gets more difficult. Hence, the tight coupling mechanism is difficult to keep high efficiency for a molecular transducer.

On the other hand, in the present loose coupling mechanism for the energy transduction, there exists chaotic motion. Then, the efficiency could be much lower than the maximum efficiency achieved by tuning initial condition in the tight coupling. However, the probability to achieve ‘good’ configurations starting from any initial conditions or from a variety of boundary conditions remains to be large. Hence energy transduction in present loose coupling mechanism works over a wide range of initial conditions, and therefore is robust.

## VI Characteristics of Energy Flow leading to Directional Motion

In spite of chaotic motion mentioned above, the output should be directional, in order to extract some work from the system. Here, recall that there is some time lag between the injection of energy and the emergence of the step motion, during which chaotic motion is maintained. Even after the stronger chaos appears with the flow of energy from the pendulum, there is a certain time lag before the emergence of the step motion (see Fig.6). During such time lag, the energy may be sustained or may be used degree by degree to bring about the step motion. In this section, we examine the energy flow in detail during this ‘active state’ preceding to the directional motion. The duration of active state is suggested to be important to realize the directional motion under large fluctuations of the configurations of molecules.

Note that the motion of the head to the neighboring lattice site is possible only if some condition for the configuration among the pendulum, head and the adjacent lattice sites is satisfied. In Fig.2, we have plotted the time course of energy flow to the chain from the head-pendulum part. Here, the time series for the kinetic energy of each lattice site is plotted together with that of the head. As shown, the energy flow from the head to farther lattice site is suppressed, before the step motion occurs. This restriction in the energy flow is generally observed whenever the directional motion is observed. On the other hand, there is a continuous gradual flow when the step motion is not possible.

To study the relevance of this restriction in the energy dissipation, we have computed energy storage by measuring the following two quantities;

$`W_{in}(t)`$ $``$ $`{\displaystyle \underset{j=i_01}{\overset{i_0+1}{}}}{\displaystyle _0^t}\dot{x}_j{\displaystyle \frac{V_\theta }{x_j}}𝑑t,`$ (5)
$`W_{out}(t)`$ $``$ $`{\displaystyle _0^t}\dot{x}_{i_02}\left\{{\displaystyle \frac{V}{x_{i_02}}}{\displaystyle \frac{K_c}{2}}(x_{i_01}x_{i_02}L)\right\}𝑑t`$ (7)
$`{\displaystyle _0^t}\dot{x}_{i_0+2}\left\{{\displaystyle \frac{V}{x_{i_0+2}}}{\displaystyle \frac{K_c}{2}}(x_{i_0+1}x_{i_0+2}L)\right\}𝑑t.`$
$`W_{in}(t)`$ is the work from the pendulum to the head and the lattice, up to time $`t`$ since the energy injection, and $`W_{out}(t)`$ is the work from the three lattice sites adjacent to the head, to farther lattice sites, again up to the time $`t`$. Here we call the head, pendulum, and the adjacent three lattice sites as ”active part”, and other lattice sites of the chain as ”residual chain part”. By adopting these terms, $`W_{out}(t)`$ is the energy dissipated to the residual chain part from the active part. In the present set of simulations for the measurement of $`W_{in}`$ and $`W_{out}`$ ( and only in this set of simulations), the contact with the heat bath to the active part is eliminated, to distinguish clearly the stored energy from the wasted work to other lattice sites.

A time course of $`W_{in}`$ and $`W_{out}`$ is shown in Figs.8. In a certain range of parameters where the step motion frequently occurs, dissipation from the active part is delayed for some time after the energy flow from the pendulum into the system (Fig.8(a)). During the interval for this delay, the energy is sustained in the active part (i.e., at the head and the adjacent lattice sites) for a while. As is shown in the left figure of Fig.8(a), the sustained energy is dissipated to the residual chain, together with the event of the step motion. (Note that steep changes in $`W_{in}`$ and $`W_{out}`$ after the step motion is an artifact since the computation here is based on the active part consisting of the head that is adjacent with the site $`i_0`$, and loses its meaning after the head steps from the $`i_0`$th site.) For the parameters with much smaller directional motion, the increase of $`W_{in}`$ and $`W_{out}`$ begins simultaneously with the energy flow from the pendulum to the system (Fig.8(b)).

The mean profile of the energy flow from the active part is given by the time course of $`dW_{out}/dt`$, where the temporal ensemble average $``$ is computed by setting the origin of $`time`$ as the time of the energy injection, since, for the parameters adopted here, the energy flows rapidly (in less than $`10`$ time units) from the pendulum to the other degrees of freedom. The rate $`dW_{in}/dt`$ takes a large value only in the initial stage ($`time<10`$) of the relaxation process. The ensemble average at time $`t`$ is computed only over the samples in which the head stays at the site $`i_0`$ up to time $`t`$. Accordingly, the number of samples decreases with the increase of the time.

As is seen in Fig.9(a), the out-flow rate $`dW_{out}/dt`$ of the energy is kept small at the early stage (up to $`time100`$) of the relaxation, for $`K_c=0.5`$. This is in strong contrast with the temporal profile for larger values of $`K_c`$ (e.g. $`=1.3`$). There, the in-flow rate $`dW_{in}/dt`$ to the active part is approximately $`0.01`$ in maximum for both the two cases. Then, the active part in the flexible system does not respond instantaneously to the sudden flow of energy, and therefore the energy is sustained in the active part without dissipation into the residual chain for a while. This sustained energy up to $`t`$ is estimated by $`W_{in}(t)W_{out}(t)`$. The ensemble average of the time course of this energy is given in Fig.9(b). The energy is sustained up to the time $`t100`$ when $`K_c=0.5`$ and dissipated more quickly when $`K_c=1.3`$.

If the injected energy were dissipated completely into a large number of degrees of freedom, the conversion of energy into mechanical work would suffer a rather large loss. Contrastingly, the motion generated by our mechanism remains confined within just a few degrees of freedom where some correlation among degrees of freedom is maintained. Supplied energy to the pendulum is not diffused as heat. As a result, the conversion is efficient. The restriction of energy transfer to the farther lattice sites comes from the correlated motion among the pendulum, head, and the lattice sites .

The effective rate for the sustained energy in the active part is estimated by $`W_{in}(t_w)W_{out}(t_w)/W_{in}(t_w)`$ by suitably choosing $`t_w`$. (As $`t_w`$ is increased to infinity, these quantities go to zero, of course). By choosing the time $`t_w`$ as an average waiting time for the step motion, one can estimate the sustained energy within the time scale for the step motion. As shown in Fig.10 approximately $`80\%`$ of the energy flow is sustained in the active part for $`K_c0.5`$, where the average steps in directional motion becomes maximum. The mean waiting time $`t_w`$ is almost $`100`$ for parameter values adapted in Fig.10, which is much shorter than the time for the complete dissipation of the energy.

Now it is shown that the combination of chaotic motion (for the loss of details of input) and the restriction of dissipation (that leads to directional motion) gives a basis of our mechanism how the energy is transduced in a microscopic system under large fluctuations of their configurations. Due to the chaotic motion, if the head shifts to the next site by crossing the energy barrier or not is probabilistic. Here, the active part can sustain the excess energy effectively. Therefore, the head comes close to the barrier several times even if it fails in the first trial. Hence the motion to cross the barrier can be robust and adaptive.

Revisiting Fig.2, it is noted that the amount of the stored energy $`E_\theta `$ in the pendulum decreases step by step, corresponding to the step motion of the head. The energy once extracted from the pendulum is partly restored after the completion of the step motion. The fact that the energy is restored to the pendulum allows the delayed use for the multiple step motion. Thus, in our mechanism the energy is used step by step.

To confirm the delayed use in multiple steps, we have measured $`W_{in}`$, by restricting only to such case that multiple steps are observed (Fig.11). Here we have computed $`W_{in}`$ up to the time of the first, second, and third steps respectively. These measurements show that the extracted energy $`W_{in}`$ from the pendulum up to the $`n`$-step motion increases with $`n`$. These results imply that the multiple step motion is supported by the delayed step-by-step use of the stored energy. This is why output motion with more and more steps is possible as the injected energy is increased.

In Fig.11, one may note that the amount of the available energy in the pendulum for the next steps is smaller than that for the previous step. Indeed, the extracted energy in Fig.11 is estimated to take a rather large value due to the timing to measure it, i.e., $`W_{in}`$ is measured at the timing when the head crosses the energy barrier to move to the next site. Then, the energy stored in the pendulum is expected to be larger than expected. This is because some part of the energy is restored to the pendulum when the head approaches the bottom of the potential for the next site, as is seen in several rises of $`E_\theta `$ after the first step motion shown in Fig.2.

## VII Discussion

In the paper, we have proposed a concept of autonomous energy transducer at a molecular scale. In this transducer, we have set the constraints (i)-(iii) in Section 2: Output is produced, even though input energy is of the same order of the thermal energy, and even though the magnitude or the timing of input is not specified, nor any control after the input is introduced. This requisite is postulated so that the function of the autonomous energy transducers is robust at a molecular scale.

We have shown a possible scheme that satisfies this requisite, by providing a simple example, with a dynamical system of several degrees of freedom. In the model, the energy transducer is robust by the following two properties: The first one is chaotic dynamics. The chaotic dynamics amplifies difference of initial conditions by each input and thus detailed condition on the input event is lost. Chaotic dynamics also brings mixing of orbits, and each orbit from different initial condition has a chance to visit a certain part of phase space that is necessary to generate directional motion. Hence, an output motion can be obtained on the average, independently of details of the input.

The second property is energy “storage” at “active part” in the system. The active part consists of several degrees of freedom, and is self-organized by absorbing energy from the input degrees of freedom. This active part is sustained over rather long time interval, with very weakly dissipating the energy. During this long time span, each orbit has more chances to visit the part of phase space inducing directional motion. Then, the ratio of the events with the directional output motion is increased.

In addition, a “storage unit”, like a pendulum, enhances the potentiality of autonomous regulation of the system. The dynamic linking between the storage and the active parts makes restorage of the extracted energy possible, which makes the delayed use for multiple steps easier. Such regulation would not be possible if the energy were simply provided to the active part, as assumed for a single conformational change per one chemical process.

As a consequence of our scheme of energy transduction, the output has large fluctuations by each event. There exists variance of the same order with the mean of the output. Hence input-output relationship is not one-to-one, but the output from the same input is distributed. Now, a mechanism is provided for a molecular machine with loose coupling between the input and the output, as was originally proposed by Oosawa.

For the autonomous energy transducer, the output has to be extracted without specifying input conditions, even under the fluctuations of the same order of the input. Hence, in order for the transducer to work at a molecular scale, several degrees of freedom are necessary to exhibit flexible change. The input energy is first stored into some of these degrees of freedom, and then is used step by step. Accordingly, the output is not directly coupled to the input. In contrast with a tight-coupling machine consisting of ”rigid gears”, the machine with this loose coupling is flexible.

After the output is extracted from the transducer, a stationary final state is restored, that is identical with the initial state before the input, except the change for the output (i.e., translational motion). In the model we presented, such stationary state is represented by a combination of stable two-body motions, that are almost decoupled each other. During the transduction process, however, several degrees of freedom in the internal states are coupled. This mode coupling leads to flexible dynamics on the one hand, but also brings about the directional motion, on the other hand. Relevance of the change of effective degrees of freedom to energy transduction has also been discussed in relationship with chaotic itinerancy .

Although we presented our scheme by a specific model with a head and a rail, the mechanism proposed here is expected to be rather general. What we have required is just a dynamical system with several degrees of freedom, including those for input and output. There, the system changes effective degrees of freedom through the internal dynamics in the system. We note here that each degree of freedom we used in the model does not necessarily have to correspond to an atomic motion. Rather, it may represent a collective variable consisting of a large number of atoms. It should be also noted that the present scheme for the energy transduction is not necessarily restricted to a Langevin system with weak damping, but it is hoped that a model with overdamped dynamics will be constructed that realizes energy transduction with loose coupling by the present scheme. In this sense, the concept of autonomous transducer at a molecular scale, which we have proposed here, would be general and applied widely.

The importance of loose coupling mechanism lies in its adaptability to different conditions. Although the motion is not precisely optimized with regards to the efficiency of the transduction, instead, the present transducer works even if external conditions are changed.

As an extreme case, consider such external condition that makes directional motion much harder. For example, by adding a load to the transducer, it is expected that the output motion is suppressed. Then the question is raised if the fluctuation in the output may also be decreased. If the answer is in the affirmative, the coupling between input and output will be more tight, as added load is increased. In other words, a transducer with loose coupling can show tight coupling depending on external condition. On the other hand, a system designed to have a tight-coupling cannot work as a machine with a loose-coupling, under the change of external condition. It just does not work when the condition is changed.

To check this plasticity of energy transducer with loose-coupling, we have carried out a numerical experiment of our model by adding a load. In order to study the effect of a load, we make a slight modification of the present model; We reverse a condition for fixing the total system, i.e., the head is fixed (e.g., by attached to a glass plate), and the rail (chain) is set free, and its center of mass can move. Due to the friction from the heat bath, the rail itself imposes a certain amount of load to the head. The amount of the load increases with the length of the rail given by $`N_c`$.

In Fig.12, the step distributions are displayed for two different lengths of chain, ($`N_c=5`$ and $`200`$). When the length of the chain is short, the mean step size is multiple that is distributed rather broadly by samples, similarly to the original system (refer Fig.7). On the other hand, for a sufficiently long chain, the mean is still kept positive but becomes smaller in amount. Moreover, the distribution gets much sharper. The fluctuation of the output is very small, i,e, the output is almost constant, as is expected in the tight coupling mechanism. Such trend is also observed in the experiments of molecular motors .

## VIII Possible Application of Autonomous Energy Transducer

A good candidate of the autonomous energy transducer is protein. In applying the present idea to proteins, it is interesting to recall that several reports show that protein is not rigid, but flexible. Actin, a rail protein, is known to be not so stiff, as assumed in the chain of our model. Some proteins are known to take several forms and spontaneously switch over these forms under thermal energy.

Recently, there are several reports on single-molecule protein dynamics, by adopting fluorescent techniques. Some enzymes show long-term relaxation of a time scale of m-sec to seconds. In particular, oscillatory dynamics with order of seconds are observed in the conformational change of single enzyme molecules, which suggest cyclic dynamics over transient states. Anomalously long-term energy storage is also observed by Ishijima et al . There the chemical energy from ATP hydrolysis is stored in the motor for 0.1-1sec and then the energy is used for the directional motion gradually. The energy is used step-by-step for multiple steps . Oscillatory dynamics are also suggested in Dynein motor . Summing up, several experiments support long-term dynamics in the relaxation of protein after excitations, as assumed in the dynamics of our model.

If the output is simply provided by chemical reaction, for instance by hydrolysis of ATP, such dynamical behavior would not be possible. For the dynamics, the proteins in concern should possess a “storage unit” between the domain for chemical reaction and the active part. In the present model, the storage unit is a simple pendulum and is not directly contacted with the heat bath. This is just an idealization, and the storage unit need not be so simple. Also, the obtained results would be valid qualitatively even though a weak contact with the heat bath is added.

From our result, it is postulated that the storage unit should interact with the active part bi-directionally, and the dissipation of energy from the unit is rather weak. Then, multiple step-by-step motion with restorage of the energy is possible. We hope that the storage unit that satisfies the above requisite is discovered in motor proteins.

Although the purpose of our paper is to propose a new general mechanism for loose coupling for energy conversion, application of the idea to a molecular motor may look straightforward. In doing so, however, one might claim that one should construct a model with overdamped dynamics, rather than a model with inertia. Although it should be important to construct an overdamped-dynamics version of the present mechanism, still the present model may be relevant to some molecular motors, for the following reasons.

First, it is not completely sure if the model for biomotor should be definitely overdamped. Proteins there exhibit fluctuations of very slow time scales, while the dynamics to realize directional motion is considerably rapid. The time scale for the stepping motion is suggested to be faster than microsecond order  and could be of nanosecond order, while the waiting time before stepping is milliseconds. In equilibrium, proteins show slow thermal fluctuations of the scale of millisecond or so . Thus the time scale of the relaxation of protein might be slower than, or of the same order as the stepping motion, in contrast with the standard estimate. It is still open if the stepping of the protein motor should be definitely treated as a overdamped motion.

It should also be noted that in the usual estimate for the damping, one uses the viscosity of macroscopic water, whereas if the water around a protein can be treated just as macroscopic water or not is an open question .

Second, as already mentioned, each variable in our model does not necessarily correspond to the atomic configuration of a protein. Rather it can correspond to some mode at a mesoscopic scale, and the inertia term in our model does not directly correspond to that of an atomic motion. When the head swings along the chain, microscopic configurations do not necessarily return to the original, but may change along the relaxation. As long as some collective mode swings, the present mechanism works.

A protein includes a large number of atoms. To understand energy transduction with such large molecule, it is important to use a reduced model with a smaller number of degrees of freedom representing a collective mode consisting of a large number of atoms.

As a theory for loose coupling, thermal ratchet models of various forms is proposed, and they have captured some features of energy conversion. However, some problems remain. First, to attain a reasonable efficiency, one needs to assume a rather tuned timing for the change of switching of potential or the time scale of colored noise. When tuned, the energy conversion often becomes one-to-one. Second, such external switching or a specific form of noise is given externally, and is not given within the theory in a self-contained form. On the other hand, in a Feynmann-type ratchet, use of different temperatures is assumed, but the temperature localized at one part of a molecule is not necessarily well defined. Some exogenous process is required. Third, the ratchet model adopts statistical description, and a single event of dynamics is not pursued, while in the recent single molecule experiments, each event of the conversion from chemical energy of ATP to mechanical work is pursued.

Whether the concept of ’autonomous energy transducer’ is valid for biomotor or not has to be judged in experiments. We discuss here possible predictions that can be tested experimentally in molecular motors.

* Mean step size per one ATP hydrolysis is modified by the change of stiffness for the rail filament. (See Fig.5 for the stiffness dependence of the mean step size). Although it might be subtle which type of stiffness should be measured to compare with our prediction, it would be important to measure stiffness dependence in any form, as a first step.
* Mean step size is modified also by the change of the amount of energy from ATP hydrolysis. If the amount of input energy is modified, the step size is predicted to change, roughly in proportion to it. Although artificial modification of the available energy by ATP is not easy experimentally, such modification of $`E_0`$, might be already realized by the long term storage of energy in molecules. If the molecule can store the energy over several times of ATP hydrolysis, the amount of available energy will simply increase. Then, the step size realized in a single sequence would be observed to be significantly large, compared with the usual step size of molecular motors. In such cases, one-to-one correspondence would be lost between ATP hydrolysis and one sequence of stepping.
* Molecular motors can step occasionally even without activation by ATP hydrolysis. A rail element which has experienced stepping of motors accompanied by ATP hydrolysis could possess long term memory and store a part of energy transduced from the motors during these events. Then, a motor attached to this rail protein would later receive this stored energy and use it for the stepping motion. Similarly, stepping events without ATP may occur due to the energy propagated from the neighboring active motors with sufficient energy from ATP hydrolysis.
* More definite proof for the relevance of autonomous dynamics in the stepping motion could be provided by direct observation of proteins dynamics. Detection of oscillation in configuration over a few periods preceding to the stepping motion will confirm the validity of the autonomous dynamics to biomotors. Measurement of conformational fluctuations with using FRETs or with using optical traps might be effective to explore this point, although for complete test of the theory, rather fine time precision, say nano or micro seconds, might be required.

To close the paper, we note that design of a nano-scale machine will be interesting that works under the present mechanism. Since our scheme is general and the model is simple, there can be several ways to realize the required situation. With the recent advances in nano-technology, the design of autonomous energy transducer at a molecular scale may be realized in near future.

Acknowledgment

We are grateful for stimulating discussion with F.Oosawa, T.Yomo, M.Peyrard and K.Kitamura. We also thank M. Peyrard for critical reading of this manuscript. The present work is supported by Grants-in-Aids for Scientific Research from Japan Society for the Promotion of Science and from the Ministry of Education, Science and Culture of Japan.


text:
# Carbon antisite clusters in SiC: a possible pathway to the D_"II" center

## I Introduction

Many unique features of SiC, such as the wide band gap and the very high electrical and thermal stability, render this semiconductor especially important for high power, high frequency and high temperature applications. As for any semiconductor material, the identification and understanding of structural defects and impurities is the key to the technological control of SiC. This is especially important in connection with the ion-implantation, which is inevitably accompanied by the generation of intrinsic defects. Some of these defects are thermally very stable. Two important examples are the photoluminescence (PL) centers D$`_\text{I}`$ (Ref. Patrick and Choyke, 1972) and D$`_{\text{II}}`$ (Ref. Patrick and Choyke, 1973). Although the centers have been extensively studied over the past 30 years,Patrick and Choyke (1972, 1973); Egilsson et al. (1999); Freitas et al. (1987); Sridhara et al. (1998, 2001); Carlsson et al. (2003) their microscopic origin remained unclear. In this work we study carbon antisite clusters and their possible relevance to D$`_{\text{II}}`$-type centers.

D$`_{\text{II}}`$ centers are always present in ion-implanted and annealed SiC samples, regardless of the implanted species and the polytype. They were also found (in low concentration) in as-grown material.Freitas et al. (1987) The abundance of D$`_{\text{II}}`$ centers increases significantly at high annealing temperatures (above 1300℃) and the center is thermally stable up to 1700℃.Freitas et al. (1987) These properties imply that the D$`_{\text{II}}`$ centers are intrinsic defect complexes. The D$`_{\text{II}}`$ luminescence originates from the recombination of an exciton bound to the defect. The spectra are usually identified by the main zero-phonon line. The important specific feature of D$`_{\text{II}}`$ are characteristic phonon replicas above the SiC phonon spectrum. In the original experiments in 3C-SiC the strongest five replicas were highlighted.Patrick and Choyke (1973) In later experimentsSridhara et al. (1998); Carlsson et al. (2003) in 4H and 6H-SiC additional localized vibrational modes (LVMs) were observed (more than 12 LVMs have been counted in 4H-SiCCarlsson et al. (2003)). Many of these observed modes are polytype-independent. A comparison of the spectra from Ref. Sridhara et al., 1998 and the earlier Ref. Patrick and Choyke, 1973 reveals that more than the five highlighted LVMs may also be contained in the 3C spectrum. Their spectral density resembles the phonon density of states in diamond, indicating that a carbon-dominated defect is responsible for the observed spectrum.

In this work we consider a number of carbon-related defects in SiC: the carbon di-interstitial, which was the first suggested model for the D$`_{\text{II}}`$ center,Patrick and Choyke (1973) the carbon antisite, the carbon split-interstitials and several carbon clusters. Comparing the calculated LVMs with the experimental D$`_{\text{II}}`$ spectra we are able to unequivocally rule out all candidates except the carbon antisite clusters. The latter are the only defects that provide a rich LVM spectrum similar to that of the D$`_{\text{II}}`$ center. The core structure of all these clusters is the dicarbon antisite (C2)$`_{\text{Si}}`$ (cf. Fig. 1 left) which acts as an aggegation center for further carbon atoms. In this paper we will first present our results for the dicarbon antisite and then consider the formation of more complex defect clusters. We calculate the LVMs associated with these defects and compare them with the D$`_{\text{II}}`$ spectrum. Although theoretically all these defects can acquire different charge states, we focus on the neutral state as the most relevant for PL. For excitons bound to a charged defect more non-radiative recombination channels are open, and they are unlikely to be seen in PL-experiments.Stoneham (1975) We also find that the formation kinetics of the carbon clusters agrees with the known features of the D$`_{\text{II}}`$ center. These results indicate that although neither of the considered defects can alone explain all features of the D$`_{\text{II}}`$ center, the carbon antisite clusters may serve as a core structure of D$`_{\text{II}}`$-type defects.

## II Method

We employ an *ab initio* density functional theory (DFT) approach as implemented in the software package FHI96SPIN.Bockstedte et al. (1997) Smooth norm-conserving pseudopotentials of the Troullier-Martin typeFuchs and Scheffler (1999) and a plane-wave basis-set with a cut-off energy of 30 Ry are used. The exchange-correlation potential is approximated within the LSDA in the parametrization of Perdew and Zunger.Perdew and Zunger (1981) In order to reduce the artificial defect-defect interaction, large supercells with 216 lattice sites for 3C-SiC and 128 sites for 4H-SiC are used for all calculations of the defect energetics. For the 216-site cell, the Brillouin zone is sampled by the $`\mathrm{\Gamma }`$ point, whereas for the hexagonal 128-site cell a special Monkhorst and PackMonkhorst and Pack (1976) $`2\times 2\times 2`$ mesh is used. The values of the formation energies are given assuming silicon-rich conditions. Under carbon rich conditions they should be reduced by roughly $`3\mathrm{\Delta }H_\text{f}=1.74`$ eV, where $`\mathrm{\Delta }H_\text{f}`$ is the heat of formation of SiC. The vibrational properties are calculated using the *frozen phonon* method. The dynamical matrix

$$\mathrm{\Phi }_{ij}=\frac{1}{\sqrt{m_im_j}}\frac{^2E}{Q_iQ_j}=\frac{1}{\sqrt{m_im_j}}\frac{F_i}{Q_j}$$

is calculated by applying displacements $`Q_j`$ from the equilibrium configuration and evaluating the force derivatives $`F_i/Q_j`$. To obtain noticeable energy changes, the displacements $`Q_j`$ ought to be much larger than realistic phonon amplitudes. Therefore it is important to eliminate the anharmonic contribution when evaluating $`\mathrm{\Phi }_{ij}`$. We do this by applying three different values of the displacement and extracting the linear force constant from a polynomial expansion of the forces. The forces have been converged to a relative accuracy of $`10^4`$. For 3C-SiC, the dynamical matrix is evaluated for a supercell with 64 or 216 sites. The experimental bulk phonon modes are reproduced within 5 percent by this method. Yet the polarization splitting of the TO and the LO modes at the $`\mathrm{\Gamma }`$ point is missing, since the macroscopic polarization of the crystal is incompatible with the periodic boundary conditions imposed on the supercell. We obtain a triple degenerate mode at 115.3 meV, in place of the TO mode of 98.7 meV and the LO mode of 120 meV.Widulle et al. (1999) For 4H-SiC and for large clusters in 3C-SiC, the computational cost of evaluating the dynamical matrix of the whole 128-site or 216-site supercell is prohibitively large. In these cases we constrain the LVM calculation to the defect molecule (i.e. the dumbbell and its nearest neighbors) embedded in a supercell. We have examined the error resulting from the defect molecule approximation and from the uncertainty in the lattice constant. The results of these tests for a dicarbon antisite in 3C-SiC (see below) are shown in Fig. 2. The lines represent the frequencies of the LVMs versus the defect molecule size calculated at the theoretical (LSDA) lattice constant. The frequencies shift within the grey area when the lattice constant is increased towards the experimental value of 8.239 Bohr. The black bars are the same frequency ranges, but in this case both the electronic structure calculation and the LVM calculation have been performed in a 64 sites cell. We have also verified that under the pressure exerted by the periodic array of defects onto the lattice the lattice constant lies between the indicated limits. The LVM’s variation due to the uncertainty in the lattice constant is thus less than 10 meV. It is also apparent that the defect molecule approximation affects the least localized modes 1 and (to a much lesser extend) 2. The largest effect on the accuracy has the inclusion of the full supercell, which shifts the mode 1 by 14 meV in the 216 sites cell. Compared to the full 64 sites cell, the mode 1 is only 2 meV higher.

## III The dicarbon antisite

Let us now consider one by one the candidates for the D$`_{\text{II}}`$ center as listed above. First, we find that the originally suggested carbon di-interstitial has a formation energy of about 12 eV. This is extremely high even for SiC and practically rules out this model. A similar defect is the carbon split-interstitial C$`_{\text{sp}100}`$, i.e. a pair of carbon atoms sharing the same (carbon) site and pointing in $`100`$ direction (in 3C-SiC). This defect also possesses a short carbon-carbon bond and it is energetically the most favorable carbon interstitial.Mattausch et al. (2001a) The second lowest carbon interstitial is the carbon-silicon split-interstitial C$`_{\text{spSi}100}`$ (a carbon and a silicon atom on a silicon site). Here the additional carbon atom is also connected to its nearest neighbors by short carbon-carbon bonds. However, according to our calculations both defects are very mobile,Bockstedte et al. (2003) which is in conflict with the thermal stability of the D$`_{\text{II}}`$ center. In addition, we find the LVMs of both defects are not compatible with the D$`_{\text{II}}`$ spectrum. The neutral carbon split-interstitial C$`_{\text{sp}100}`$ in the tilted structureBockstedte et al. (2003) has two LVMs at 120.4 and 188.7 meV and LVMs in the phonon gap. The carbon-silicon split-interstitial C$`{}_{}{}^{2+}{}_{\text{spSi}100}{}^{}`$ possesses a richer vibrational spectrum but still with only three non-degenerate LVMs above the bulk spectrum. This also does not match the number of the characteristic D$`_{\text{II}}`$ lines. Another possible candidate could be the carbon antisite C$`_{\text{Si}}`$. Following the carbon vacancy V$`_\text{C}`$ it is the second most abundant defect in SiC. However, it shows only vibrational resonances in the bulk phonon spectrum and has no true LVMs. Also small clusters of carbon interstitials do not provide LVMs compatible with the D$`_{\text{II}}`$ spectrum.Mattausch et al. (2003)

The simplest candidate that possesses the main features of the D$`_{\text{II}}`$ center is the dicarbon antisite (C2)$`_{\text{Si}}`$ that we have already mentioned above in connection with Fig. 2. The dumbbell of two carbon atoms sharing a silicon site (highlighted atoms in Fig. 1 left) is surrounded by a tetrahedron of four carbon atoms. In 3C-SiC the dumbbell is pointing in $`100`$ direction whereas in 4H-SiC it is oriented either along the $`10\overline{1}1`$ or the $`\overline{1}011`$ direction. Depending on the Fermi level, the defect may acquire the charge states from $`2^+`$ through $`2^{}`$.

We find that the dicarbon antisite has a relatively low formation energy. For the neutral charge state the values are 6.7 eV and 7.8 eV in 3C and 4H polytypes, respectively. In 4H-SiC two possible substitutional configurations exist, i.e. hexagonal (*h*) and cubic (*k*), but the formation energy at these two sites only differs by 0.1 eV. These values are comparable to the formation energies of C$`_{\text{sp}100}`$ and C$`_{\text{spSi}100}`$. Still, the calculated formation energy indicates that its equilibrium abundance is too low for experimental observation. Hence kinetic effects should play an important role in the dicarbon antisite formation. On the other hand, once created, the defect has a very high stability, as reflected by the high binding energy. Before we turn to kinetic aspects, we discuss the electronic properties and the LVMs of the dicarbon antisite.

In 3C-SiC, the dicarbon antisite possesses two degenerate defect orbitals with energy levels in the band gap. In the neutral charge state these orbitals are occupied by two electrons, which can form either a spin 0 (low-spin) or a spin 1 (high-spin) configuration. According to the calculation, the high-spin state is about 80 meV lower in energy than the low-spin state. This energy difference is, however, too small to uniquely determine the ground state spin configuration. For this reason we consider both options. The defect wave function consists of bonds between the three upper carbon atoms and the three lower carbon atoms of the defect tetrahedron (cf. Fig. 1 left). In the low-spin configuration, one of these bonds is fully occupied, resulting in a splitting of the defect levels and a Jahn-Teller distortion, with the upper or the lower pair of the tetrahedron being stretched. This reduces the original $`D_{2d}`$ symmetry to $`C_{2v}`$. In the high-spin configuration, both bonds are equally occupied, retaining the degeneracy and the $`D_{2d}`$ symmetry. As we discuss below, this geometrical difference is reflected in the vibrational modes of the defect. In 4H-SiC, the high-spin state is about 120 meV lower than the low-spin state, which is—as in 3C-SiC—a very small energy difference. Due to the lower symmetry of the crystal, the symmetry group is always $`C_{1h}`$. The defect’s structure on the cubic and the hexagonal site is qualitatively similar.

As outlined above, the dicarbon antisite is the simplest defect that possesses LVMs that are similar to the phonon spectrum of the D$`_{\text{II}}`$ center. The calculated frequencies (using a 64-sites supercell) for the neutral charge state in low-spin and high-spin configurations are listed in Tab. 1 along with the experimental data for the five characteristic D$`_{\text{II}}`$ phonon modes and our calculated results for the two neighboring dicarbon antisites. We see that in the 3C polytype the low spin configuration provides an LVM pattern similar to the five phonon replicas of D$`_{\text{II}}`$. Similar results are obtained for 4H-SiC, except that the frequencies of the modes 1 to 4 are lowered due to the defect molecule approximation.

In spite of a different crystal symmetry, the vibrational patterns of the LVMs in 3C- and 4H-SiC are very similar, owing to the localized nature of the vibrational modes. The highest mode 5 is a stretching vibration, where the two atoms of the dumbbell oscillate against each other. Its calculated frequency is about 8% higher than the experimental value. However, this is the usual accuracy of these defect calculations.Capaz et al. (1998) The modes 3 and 4 represent an oscillation of one of the dumbbell atoms against its neighbors. In 3C-SiC, these modes are affected by the Jahn-Teller distortion of the low-spin configuration. As expected from the symmetry, the modes are degenerate in the high-spin and non-degenerate in the low-spin case. This effect can also be seen in Fig. 2, where the 216 sites calculation, although performed spin-unpolarized, has been performed in a high-spin like configuration. Considering the small energy difference between the high-spin and the low-spin configuration both options should be allowed. In addition to the already mentioned clustering of carbon atoms, this may be a source of further lines in the luminescence spectra. In 4H-SiC, the modes 3 and 4 are always non-degenerate, although the splitting in the high-spin configuration is smaller than in the low-spin configuration, due to the lower symmetry of the crystal field. The mode 2 describes the vibration of the dumbbell along its axis against the enclosing tetrahedron. The mode 1 is the least localized mode and therefore it is most affected by the employed supercell technique. This mode corresponds to a breathing-like vibration of the enclosing tetrahedron. While for the low-spin case the breathing mode lies above the highest bulk vibration, it is in resonance with the bulk spectrum and lies 8 meV below the calculated value of the highest bulk mode for the high-spin case in a 64 sites cell. This is in agreement with Gali *et al.* (Ref. Gali et al., 2002) who considered the high-spin state of (C2)$`_{\text{Si}}`$. In 4H-SiC, a similar behavior is not obvious from our results, as the effect may be masked by the defect molecule approximation. In addition to the high-frequency LVMs discussed above, we also find maxima of the vibrational density of states in the phonon band gap.

It is noteworthy that the results do not change significantly for different possible charge states of the defect.Mattausch et al. (2001b); Mattausch et al. (2002) We observe a slight lowering of the LVM’s frequencies with the increase of the number of localized electrons, resulting from an expansion of the defect molecule due to the occupation of the bonds. For 3C-SiC, the modes 3 and 4 are degenerate in the charge states $`2^+`$ and $`2^{}`$, since the electronic defect levels are either empty or fully occupied, and a Jahn-Teller distortion is suppressed.

## IV Clustering of carbon atoms

Usually, the D$`_{\text{II}}`$ center is observed in ion-implanted and annealed samples. It is stable up to temperatures of 1700℃. This and the previously discussed low equilibrium concentration suggest a kinetically controlled formation of this defect from non-equilibrium vacancies and interstitials. At first, a carbon antisite is formed via a recombination of a carbon interstitial with a silicon vacancy. It merges then with a mobile carbon split-interstitial, producing a dicarbon antisite with a high binding energy that ranges from 3.9 eV ((C2)$`{}_{}{}^{2+}{}_{\text{Si}}{}^{}`$) to 5.1 eV ((C2)$`{}_{}{}^{2}{}_{\text{Si}}{}^{}`$). The dicarbon antisite can serve as a condensation center for larger clusters. Such clusters grow by aggregating carbon split-interstitials around the core. Adding a further carbon atom to the dicarbon antisite releases a binding energy of about 4.8 eV. Binding one more carbon interstitial to this complex yields a very small cluster of four carbon atoms at a silicon site (cf. Fig. 1 right). According to our calculations the energy of 2.8 eV is needed to remove a carbon atom from this structure. However, for this defect we obtain LVMs well above the D$`_{\text{II}}`$ spectrum. A detailed discussion of carbon clusters will be given elsewhere.Mattausch et al. (2003)

In addition to the clusters located around a single silicon site, aggregates of several dicarbon antisites are also possible. A pair of the dicarbon antisites in 3C-SiC ((C2)$`_{\text{Si}}`$)2 is shown in Fig. 3. In 4H-SiC three different configurations, which can be compared to its 3C-SiC counterpart, are possible when the two di-antisites occupy two cubic (*kk*), two hexagonal (*hh*) or cubic-hexagonal (*kh*) sites. As an exemplary structure we have considered in detail the *kh*-configuration, which has a local geometry similar to the 3C-SiC case (Fig. 3), but the bond length between the neighboring dicarbon antisites is shorter. We found that this defect is very stable in both polytypes: in 3C-SiC it costs about 5.9 eV to remove a carbon atom, whereas in 4H-SiC this energy is 6.7 eV. The calculated LVMs of both complexes are given in Tab. 1. They have been calculated using the defect molecule approximation that allows the vibration of the two dumbbells and of the surrounding carbon neighbors within the 216 (3C-SiC) and 128 (4H-SiC) sites cell. As discussed above, the numerical results are affected by the defect molecule approximation. Alternatively, one could work with the smaller 64 sites cell which would allow to calculate the dynamical matrix of the whole cell. However, we find that in this cell the electronic structure is so strongly distorted by the proximity effects that it is unacceptable for the analysis of the LVMs. For the pair of dicarbon antisites the highest almost degenerate modes 5 and 6 involve a stretching vibration of the dumbbells, similar to the lone dicarbon antisite. Due to the interaction of the dumbbells these modes become softer and combine into anti-symmetric and symmetric vibrations. The modes 3 and 4 represent a motion of the two lower and upper carbon atoms of the dumbbells in bond direction against their surrounding neighbors, respectively. In 4H-SiC, the defect shows a similar behavior, but the energy of all modes is by about 8 meV lower than in 3C-SiC so that the lowest mode (with an energy of 114.4 meV) drops into the bulk phonon spectrum. This effect and also the larger frequency difference between the stretching modes reflects the stronger relaxation of the defect in this polytype. We have recently learned that A. Gali *et al.* calculated the LVMs of the $`kk`$-configuration, finding results similar to the 3C-SiC values, and speculated that this defect was the D$`_{\text{II}}`$ center.Gali et al. (2003) As the $`kk`$-configuration has a 3C-like environment (in contrast to the $`kh`$-configuration), the similarity to 3C-SiC is quite natural. However, the conjecture that the ((C2)$`_{\text{Si}}`$)2 model completely describes the D$`_{\text{II}}`$ center does not seem to be supported by experimental data. Recent experiments have investigated the multitude of zero phonon lines (ZPLs)Sridhara et al. (2001) as well as the phonon replicasCarlsson et al. (2003) of the D$`_{\text{II}}`$ center. In 6H-SiC, Sridhara *et al.*Sridhara et al. (2001) concluded from temperature and stress-dependence experiments that the four ZPLs of D$`_{\text{II}}`$ stem from the different excitonic states of the bound exciton (one ground state and three excited states), but not from different binding centers of the same exciton. Regarding the phonon replicas in 4H-SiC more than 12 LVMs and a single ZPL have been reported by Carlsson *et al.*,Carlsson et al. (2003) which is consistent with the earlier results for 4H and 6H-SiC of Sridhara *et al*.Sridhara et al. (1998) In addition, as we noted above, a comparison of the spectra in 3C,Patrick and Choyke (1973) 4H and 6H-SiCSridhara et al. (1998) indicates that more than five LVMs may be present also in 3C-SiC. The apparent large number of phonon replicas in 4H and 6H-SiC makes it tempting to assume that they originate from different geometrical configurations of the same simple defect. For example, for the ((C2)$`_{\text{Si}}`$)2 in 4H-SiC this would be $`kk`$, $`hh`$ and $`kh`$. This assumption has, however, the following implications. Firstly, it follows that a larger number of LVMs should be present in 6H-SiC than in 4H or 3C-SiC, since less or no inequivalent sites are available in the two latter polytypes. Secondly, the different configurations should give rise to a multitude of exciton ground-state ZPLs. According to the conclusion of Sridhara *et al.*Sridhara et al. (2001) and the rich vibrational structure contained in all polytypes (cf. discussion above) this is not evident from the experiments. This applies, in particular, to the ((C2)$`_{\text{Si}}`$)2 defect. Furthermore, the LVM pattern of this defect is, in fact, not compatible with the D$`_{\text{II}}`$ spectrum. The small splitting of the two highest LVMs (see Tab. 1), which persists through different polytypes and configurational models, is not observed in the D$`_{\text{II}}`$ spectrum which has nearly equidistant lines. Thus the ((C2)$`_{\text{Si}}`$)2 model alone cannot explain all features of the D$`_{\text{II}}`$ center. Larger carbon clusters are likely to be involved. For such clusters the differences between the polytypes lose their importance. This may be responsible for the polytype-independent LVMs of the D$`_{\text{II}}`$ centers. Further experimental data are necessary to clarify the situation, especially for the 3C polytype, which lacks inequivalent lattice sites.

## V Conclusion

In conclusion, we have analyzed the structural, electronic and vibrational properties of the dicarbon antisite complex. The defect possesses a characteristic vibrational pattern in the frequency range of the D$`_{\text{II}}`$ spectrum, but cannot explain all the phonon replicas observed. Yet, we find that the dicarbon antisite can act as an aggregation center for carbon clusters that allow for a much richer phonon spectrum. All the defect complexes presented above possess very high dissociation energies, which guarantee their high thermal stability. Due to the high formation energies, these defects can exist in as-grown material only in very low concentrations but can be kinetically created during the annealing. Yet, none of the considered defect complexes can explain all features of the D$`_{\text{II}}`$ center. Most probably the D$`_{\text{II}}`$-type centers are related to larger carbon clusters with the dicarbon antisite playing the role of an elementary building block of the aggregates.

###### Acknowledgements.

We acknowledge fruitful discussions with W.J. Choyke. This work has been supported by the Deutsche Forschungsgemeinschaft within the SiC Research Group.