text:
# Operating a phase-locked loop controlling a high-Q tuning fork sensor for scanning force microscopy

## I introduction

Since its invention by Binnig, Quate and Gerber , the scanning force microscope (SFM) has become a standard tool for the investigation of conducting and insulating surfaces on the atomic scale. Various imaging modes exist that can be classified to be either contact modes or dynamic modes . The latter are characterized by the fact that the tip oscillates above the surface and the parameters of the oscillation are influenced by the tip-sample interaction. True atomic resolution has been achieved in this imaging mode . Usually the cantilever oscillation is measured by optical means. Tuning fork sensors offer the possibility of non-optical detection of the cantilever oscillation via the piezoelectric effect . It was demonstrated that atomic resolution is possible with these unconventional and very stiff sensors . The fundamental limits to force detection with quartz tuning forks were disussed by Grober and co-workers in Ref. .

We employ the tuning fork sensors in a cryo-SFM setup in which the home-built microscope is operated in a 3He-cryostat at a temperature of 300 mK. Tuning fork sensors can be classified by the strength of the mixing between symmetric and antisymmetric mechanical oscillation modes. In the extreme case, one of the tuning fork prongs is firmly attached to a support  and extremely strong mode mixing occurs. In the first section of this paper we present evidence, that our sensors are in the weak mode mixing limit and that the simple harmonic oscillator model is therefore appropriate. In the second section we describe our nested SFM feedback system comprising the tuning fork sensor, a phase-locked loop and a conventional $`z`$-feedback. We present an analysis of the linear response and the noise of these feedbacks and show how the optimum settings for the feedback parameters can be found. These optimum parameters are found to be independent of the tip shape and the details of the tip sample interaction. The advantages of using a feedback with phase-locked loop for the very high $`Q`$ tuning fork sensors are discussed and a comparison is made with the phase control mode.

## II Tuning fork sensors: a realistic model

In our dynamic mode SFM we employ the same type of tuning fork sensors previously discussed and calibrated in Refs. . Figure 1 shows a typical sensor. The commercially available tuning fork  is mounted on a print plate under an angle of about 10. A 15 $`\mu `$m PtIr wire is glued  to the end face of one tuning fork prong and to a copper post which serves as the electrical contact to the tip. The free end of the wire is then etched electrochemically resulting in a sharp tip which will serve as a probe for the sample surface. The preparation technique impairs the symmetry of the tuning fork and we have to care about the effects on the oscillation properties.

Figure 2 shows the current through the tuning fork near resonance measured at a temperature of 4.2 K with an excitation voltage of $`U_0=10`$ $`\mu `$Vrms applied to the tuning fork contacts. The top graph is the magnitude of the current. The experimental details of the admittance measurements will be discussed later in the paper. In the following, we start with a mechanical model of the tuning fork sensor which eventually aims at the detailed understanding of such resonance curves. In particular, we will focus on the symmetry breaking effect of the tip prepraration.

We describe the mechanical oscillator with a coupled harmonic oscillator model driven by the voltage $`U_0`$ via the piezoelectric coupling constant $`\alpha `$.

$`m_1\ddot{x}_1+\gamma _1\dot{x}_1+k_1x_1+k_cx_2+\gamma _c\dot{x}_2`$ $`=`$ $`\alpha U_0`$ (1)
$`m_2\ddot{x}_2+\gamma _2\dot{x}_2+k_2x_2+k_cx_1+\gamma _c\dot{x}_1`$ $`=`$ $`+\alpha U_0+F_{ts}(x_2).`$ (2)
The deflection of the two prongs is described by the coordinates $`x_i`$. The $`m_i`$ are the effective masses of the two prongs, the $`k_i`$ and $`\gamma _i`$ are their spring- and damping-constants, respectively. The tip is attached to prong 2 and interacts with the sample surface via the interaction force $`F_{ts}(x_2)`$. Coupling between the prongs is described by coupling constants $`k_c`$ and $`\gamma _c`$.

The Eigenmodes of this oscillator (determined with $`F_{ts}(x_2)=0`$) are given by

$`\mathrm{\Omega }_{1,2}^2=\omega _0^2\pm \omega _c^2\sqrt{1+\kappa ^2}`$

where $`\omega _0^2=(k_1/m_1+k_2/m_2)/2`$, $`\omega _c^2=k_c/\sqrt{m_1m_2}`$ and $`\kappa =(k_1/m_1k_2/m_2)/(2k_c/\sqrt{m_1m_2})`$. In the case of $`\kappa =0`$, in particular when the two prongs are identical, the two Eigenmodes are the symmetric (the $`x_i`$ oscillate with zero phase difference) and the antisymmetric (the $`x_i`$ oscillate with a phase difference of $`\pi `$) mode. The metallic tuning fork contacts are arranged such that only the antisymmetric mode is excited by a driving voltage. For $`\kappa 1`$ the two modes do not mix significantly and the symmetries of the modes are not strongly affected, while for $`\kappa 1`$ each prong has its own resonance frequency where the other is hardly excited.

In the following we show that our tuning fork sensors are in the regime of negligible mode mixing, because the introduced asymmetry is small, i.e. $`\kappa 1`$. We can rewrite the expression for $`\kappa `$ introducing $`m_1=m`$, $`m_2=m+\mathrm{\Delta }m`$, $`k_1=k`$, $`k_2=k+\mathrm{\Delta }k`$. Expanding up to first order in $`\mathrm{\Delta }m/m`$, $`\mathrm{\Delta }k/k`$ and $`k_c/k`$ we obtain

$`\kappa ={\displaystyle \frac{\mathrm{\Delta }m/m\mathrm{\Delta }k/k}{2k_c/k}}`$

It can be seen that the additional mass on one prong can be compensated by increasing the spring constant, a fact exploited already in Ref. .

We estimate the additional spring constant $`\mathrm{\Delta }k`$ of a thin wire attached to one prong. It is given by  $`\mathrm{\Delta }k=3\pi Er^4/(4L^3)`$, where $`L`$ denotes the length of the wire, $`E`$ is Young’s modulus of the material and $`r`$ is its radius. Typical wires have $`r=7.5`$ $`\mu `$m, $`L=0.5`$ mm, and $`E=200400`$ GPa leading to $`\mathrm{\Delta }k=1224`$ N/m. With $`k=14000`$ N/m this results in $`\mathrm{\Delta }k/k10^3`$.

An estimate of the additional mass added to one arm of the tuning fork takes essentially the glue into account. The idealized drop of glue has the shape of a half-sphere with radius $`r`$ and volume $`V=2\pi /3r^3`$. Using the mass density of the H20E-epoxy $`\rho =2600`$ kg/m3 and a radius of 100 $`\mu `$m we obtain a mass of 5.44 $`\mu `$g. This weight has to be related to the mass of a single tuning fork arm which is 1.1 mg. This leads to the ratio $`\mathrm{\Delta }m/m=5\times 10^3`$.

The ratio $`k_c/k`$ can be estimated from finite element simulations of tuning fork sensors to be of the order of several percent .

From these estimates it seems to be reasonable that the employed sensors are in the limit of small mode mixing. An additional observation strongly supports this point of view: only very rarely is a second resonance observed in the admittance. Since the measured current detects only the antisymmetric component of the oscillation, strong mode mixing would make the second resonance visible.

Motivated by these estimates we approximate the equations of motions by $`m_1=m_2=m`$, $`k_1=k_2=k`$ and $`\gamma _1=\gamma _2=\gamma `$ and obtain

$`m\ddot{\xi }+(\gamma \gamma _c)\dot{\xi }+(kk_c)\xi `$ $`=`$ $`2\alpha U+F_{ts}(s+\xi /2)`$ (3)
$`m\ddot{s}+(\gamma +\gamma _c)\dot{s}+(k+k_c)s`$ $`=`$ $`{\displaystyle \frac{1}{2}}F_{ts}(s+\xi /2)`$ (4)
with $`s:=(x_1+x_2)/2`$ is the centre of mass coordinate and $`\xi =x_2x_1`$ is the relative coordinate. The tip-sample interaction force couples the centre of mass motion to the relative motion. However, due to the extremely high stiffness of our tuning forks (of the order of $`10^4`$ N/m) a typical oscillation amplitude of 1 nm gives rise to a restoring force of $`10^5`$ N. This is very large compared to typical tip-sample interaction forces of the order of 1 nN. Therefore we neglect the coupling and describe the tuning fork oscillation with the approximation $`s=0`$ which implies that $`x_1=x_2`$ and $`\xi =2x_2`$. Inserting these relations into equation (3) leads to the single harmonic oscillator approximation for our tuning fork sensors, which we rewrite with $`x=x_2`$, $`(kk_c)/m=\omega _0^2`$ and $`(\gamma \gamma _c)/m=\omega _0/Q`$

$`\ddot{x}+{\displaystyle \frac{\omega _0}{Q}}\dot{x}+\omega _0^2x={\displaystyle \frac{\alpha }{m}}U+{\displaystyle \frac{1}{m}}F_{ts}(x).`$

Summarizing the discussion of the tuning fork sensors, we have found two reasons why the single harmonic oscillator description is appropriate for our sensors: first, the added mass (the glue) is small compared to the mass of a single tuning fork arm. Second, the relative mechanical coupling strength of the two prongs is strong compared to the influence of the added mass and the added spring constant. In this respect these sensors are completely equivalent to other sensors with one prong firmly attached to the support which act essentially as an extremely stiff piezoelectric cantilever. In other respects there are important differences: we find quality factors $`Q`$ of up to 250000 under the UHV-conditions occuring in our evacuated sample space at a temperature of 300 mK. These values are 1-2 orders of magnitude larger than those reported for the other type of tuning fork sensors . The implications of this will be discussed in the remainder of the paper.

## III Admittance measurement

### A Model for the admittance

The mechanical oscillation of tuning fork sensors is measured via the piezoelectric effect of the quartz crystal. The induced piezoelectric charge on one tuning fork electrode is given by $`Q_p=\alpha \xi =2\alpha x`$ and the corresponding piezoelectric current is $`I_p=2\alpha \dot{x}`$. In addition, a current of size $`I_c=i\omega CU`$ flows through the capacitance between the tuning fork contacts $`C_0`$ and the total current is therefore given by $`I(\omega )=i\omega C_0U+2i\alpha \omega x(\omega )`$. If we neglect the tip-sample interaction force we can determine $`x(\omega )`$ from the harmonic oscillator model and find

$`I(\omega )=i\omega U\left(C_0+{\displaystyle \frac{2\alpha ^2/m}{\omega _0^2\omega ^2+i\omega \omega _0/Q}}\right).`$

Resonance curves like the one shown in Fig. 2a) can be excellently fitted with this equation.

Near resonance ($`\omega \omega _0`$) the current is dominated by the piezoelectric contribution for our high $`Q`$ sensors and can be expanded to be

$$I(\omega )=\frac{2\alpha ^2QU}{m\omega _0}\left(1i\frac{2Q}{\omega _0}(\omega \omega _0)\right)+𝒪\left[(\omega \omega _0)^2\right].$$
(5)
The component of the current which is shifted by 90 with respect to the driving voltage is depicted in Fig. 2b). It is proportional to the deviation of the excitation frequency from the resonance frequency. This fact is utilized for the AFM-feedback.

### B Admittance measurement

We measure this current by using an I-U converter with a current to voltage conversion ratio of $`K=10^6`$ V/A at 32 kHz. A guard driver neutralizes the huge capacitance $`C_K=1.8`$ nF of the long coax-cable connecting the tuning fork in the cryostat and thereby increases the bandwidth of the I-U converter and reduces the output noise. It turns out that the output noise of the converter which was found in this setup to be $`n_{IUC}=800`$ nV/$`\sqrt{\text{Hz}}`$ at 32 kHz, dominates the noise in the whole AFM feedback.

Figure 3 shows the setup for the admittance measurement of the tuning fork. The fork (TF) is driven by voltage controlled oscillator (VCO) with typical excitation amplitudes between 10 $`\mu `$V and 1 mV depending on the $`Q`$-value of the tuning fork. The tuning fork current is converted into a voltage with the I-U converter (IUC). The resonance frequency of the tuning fork sensor depends on the tip sample interaction via the tip-sample separation $`\mathrm{\Delta }z`$ and the tip-sample voltage $`U_{ts}`$ (see below).

The resulting voltage is demodulated by the lock-in amplifier with sensitivity $`S`$ which determines the in-phase (X) and 90 (Y) current components corresponding to the real and imaginary components in eq. (5). The output voltages up to first order in the frequency shift are given by

$`X_{res}`$ $`=`$ $`{\displaystyle \frac{10\text{V}}{S}}K{\displaystyle \frac{\alpha ^2QU}{\pi mf_0}}`$
$`Y_{res}`$ $`=`$ $`{\displaystyle \frac{10\text{V}}{S}}K{\displaystyle \frac{\alpha ^2Q^2U}{\pi ^2mf_0^2}}\mathrm{\Delta }f=X_{res}{\displaystyle \frac{Q}{\pi f_0}}\mathrm{\Delta }f,`$
where we have introduced $`f_0=\omega _0/(2\pi )`$ and the frequency shift $`\mathrm{\Delta }f=ff_0`$. The $`Y`$-component of the signal is a linear indicator of the deviation from the resonance frequency which can be utilized for the AFM-feedback.

The spectral noise density on the output signals of the lock-in is given by $`n_Y=10\text{V}/Sn_{IUC}`$. The integrated output noise $`\delta Y_{res}`$ depends on the lock-in time constant $`\tau `$ and on the order $`n`$ of its low-pass filter via

$`\delta Y_{res}=c_n\times {\displaystyle \frac{10\text{V}}{s}}\times \sqrt{{\displaystyle \frac{n_{IUC}^2}{2\pi \tau }}}.`$

The constants $`c_n`$ are calculated to be $`c_1=1.25331`$, $`c_2=0.886227`$, $`c_3=0.767495`$, $`c_4=0.700624`$, $`c_5=0.655374`$, $`c_6=0.621742`$. In our setup we use a lock-in time constant of $`\tau =100`$ $`\mu `$s corresponding to a bandwidth of 1.6 kHz and a fourth order low-pass filter. With $`S=10`$ mV this leads to $`n_Y=800`$ $`\mu `$V/$`\sqrt{\text{Hz}}`$ and to typical values of $`\delta Y_{res}=22`$ mV corresponding to a frequency noise of 0.3 mHz.

### C Frequency shift and tip-sample interaction

The harmonic oscillator approximation for the tuning fork sensors allows us to apply the Hamilton-Yacoby perturbation result for the frequency shift

$`\mathrm{\Delta }f(\mathrm{\Delta }z,U_{ts})=`$ (7)
$`{\displaystyle \frac{f_0}{2\pi Ak}}{\displaystyle _0^{2\pi }}𝑑xF_{ts}(\mathrm{\Delta }z+A\mathrm{sin}x,U_{ts})\mathrm{sin}x.`$
$``$ $`\alpha \mathrm{\Delta }U_{ts}+\beta \mathrm{\Delta }z,`$ (8)
where $`\alpha =\mathrm{\Delta }f(0,U_0)/U_{ts}`$ and $`\beta =\mathrm{\Delta }f(0,U_{ts})/\mathrm{\Delta }z`$. The approximate expression of the frequency shift has been obtained by an expansion to first order in the tip-sample separation $`\mathrm{\Delta }z`$ and the tip-sample voltage $`U_{ts}`$ which generates electrostatic tip-sample forces. The action of these quantities on the admittance measurement are indicated in Fig. 3.

## IV Demodulation of the tuning fork sensor on resonance

### A Step Response of the Sensor

For the AFM-operation it is crucial to determine the bandwidth of the tip-sample interaction sensor, i.e., the time it takes to respond to a given step in $`\mathrm{\Delta }z`$ or $`U_{ts}`$. We first determine the step response of our sensor similar to Albrecht . Following eq. (3) we assume that the oscillator is at all times driven by the force $`2\alpha U_0\mathrm{cos}(\omega _et+\mathrm{\Phi })`$. For $`t<0`$ we excite on resonance, i.e. $`\omega _e=\omega _0`$ and assume that the oscillator performs its steady state oscillations. At $`t=0`$ we instantaneously change the resonance frequency of the oscillator to the new frequency $`\mathrm{\Omega }_0=\omega _0+\mathrm{\Delta }\omega `$ and seek solutions of eq. (3), neglecting $`F_{ts}`$, of the form

$`\xi (t)=Y(t)\mathrm{cos}(\omega _et+\mathrm{\Phi })+X(t)\mathrm{sin}(\omega _et+\mathrm{\Phi })`$

For $`t<0`$ the Amplitudes are $`X(t)=A_<`$ and $`Y(t)=0`$. For $`t>0`$ we find

$`Y(t)`$ $`=`$ $`A_>\mathrm{cos}\varphi +Be^{\omega _0t/(2Q)}\mathrm{cos}[\mathrm{\Delta }\omega t+\theta \mathrm{\Phi }]`$ (9)
$`X(t)`$ $`=`$ $`A_>\mathrm{sin}\varphi Be^{\omega _0t/(2Q)}\mathrm{sin}[\mathrm{\Delta }\omega t+\theta \mathrm{\Phi }]`$ (10)
The constants $`B`$ and $`\theta \mathrm{\Phi }`$, found from the boundary conditions at $`t=0`$, are given by

$`B^2`$ $`=`$ $`A_>^2+A_<^2+2A_>A_<\mathrm{sin}\varphi `$ (11)
$`\mathrm{cot}(\theta \mathrm{\Phi })`$ $`=`$ $`{\displaystyle \frac{A_>\mathrm{cos}\varphi }{A_<+A_>\mathrm{sin}\varphi }}`$ (12)
So far the results are generally valid and no approximations have been made. Now we consider approximate behaviour. For small frequency steps $`\mathrm{\Delta }\omega \omega _0/(2Q)`$ we find that

$`Y(t)=2A_<{\displaystyle \frac{Q\mathrm{\Delta }\omega }{\omega _0}}\left[1e^{\omega _0t/(2Q)}\mathrm{cos}\left(\mathrm{\Delta }\omega t+{\displaystyle \frac{2Q\mathrm{\Delta }\omega }{\omega _0}}\right)\right]`$

Since the long-time behaviour of the response will be governed by the exponential decay we may approximate the response by

$`Y(t)=2A_<{\displaystyle \frac{Q\mathrm{\Delta }\omega }{\omega _0}}\left[1e^{\omega _0t/(2Q)}\right],`$

i.e., a first order low-pass behaviour with bandwidth $`f_G=f_0/(2Q)`$. This gives immediately the step response of the demodulated current

$`Y_{res}(t)=X_{res}{\displaystyle \frac{Q}{\pi f_0}}\mathrm{\Delta }f\left[1e^{\pi f_0t/Q}\right],`$

or, after Fourier transformation

$`Y_{res}(\omega )=k(\omega )\mathrm{\Delta }f=X_{res}{\displaystyle \frac{Q}{\pi f_0}}\mathrm{\Delta }f{\displaystyle \frac{1}{1+i\omega /(2\pi f_G)}}.`$

### B Measurement of the Sensor Response

A measurement of this low-pass behaviour is shown in Fig. 4. The measurement was performed at a temperature of 4.2 K with the tip at a distance of 80 nm from the surface of a GaAs/AlGaAs heterostructure in which a two-dimensional electron gas (2DEG) resides 34 nm below the surface. The tuning fork was driven on resonance with an AC-voltage of 10 $`\mu `$V applied to the tuning fork contacts. The 2DEG was employed as the metallic counter electrode of the tip. A DC tip-sample voltage of $`U_{ts}=9`$ V was applied to the 2DEG while keeping the tip grounded. Following the idea of eq. (8) a low-frequency AC-voltage $`\mathrm{\Delta }U_{ts}=300`$ mV was added in order to modulate the resonance frequency of the tuning fork. Fitting the low-pass behaviour (not shown) results in the characteristic frequency of the tuning fork $`f_G=0.13`$ Hz corresponding to $`Q=116000`$. In other experiments at low temperatures we achieved $`Q`$-values of up to 250000.

In order to summarize the results of this section we want to emphasize the following points: the $`Y`$-output of the lock-in amplifier is used as the frequency detector in our setup. For small changes in resonance frequency $`\mathrm{\Delta }ff_G=f_0/(2Q)`$ the response shows low pass behaviour. Larger changes in resonance frequency lead to an oscillatory step response of the phase signal. High $`Q`$ sensors have a stronger tendency to an oscillatory response due the their higher $`f_G`$ as compared to lower $`Q`$ sensors.

## V FM-detection with a phase-locked loop

FM detection was introduced in dynamic force microscopy by Albrecht and co-workers . A theoretical description is due to Dürig and co-workers . In such setups an FM demodulator is used for measuring the oscillation frequency of the cantilever. Phase-locked loops are one among many ways of FM demodulation. They are an established method of frequency measurement . The phase-locked loop was first employed in scanning probe measurements by Dürig and co-workers . A digital version of a PLL has been developed by Ch. Loppacher and co-workers .

Due to the high $`Q`$-values of our sensors the response to changes in resonance frequency, in particular during the approach of the tip to the surface, is too slow to allow reasonable approach speeds. This is one reason why it is advantageous to employ a phase-locked loop in our setup. As we will show below, the phase-locked loop (PLL) increases the bandwidth of the response at the expense of increased noise.

### A Linear response

The PLL-setup is schematically shown in Fig. 5. The output of the interaction sensor (IS), i.e., the phase signal $`Y`$, is fed into a P-I controller (PI). Its output $`V_{\mathrm{\Delta }f}`$ is used as the input of the VCO in the interaction sensor. In the following we describe the PLL by the linear equations

$`Y(\omega )`$ $`=`$ $`k(\omega )(f_{res}(\omega )f_d(\omega ))`$ (13)
$`V_{\mathrm{\Delta }f}(\omega )`$ $`=`$ $`R(\omega )Y(\omega )`$ (14)
$`f_d(\omega )f_0`$ $`=`$ $`\eta V_{\mathrm{\Delta }f}(\omega ),`$ (15)
where $`f_{res}`$ is the resonance frequency of the tuning fork given by external parameters such as $`\mathrm{\Delta }z`$ and $`U_{ts}`$, $`f_d`$ is the frequency output of the VCO driving the tuning fork and $`f_0`$ is the center frequency of the VCO around which $`f_d`$ can be modulated. The function $`k(\omega )`$ is the response of the interaction sensor shown in Fig. 3 including the lock-in demodulator with its frequency response defined by the lock-in time constant and filter characteristics, $`R(\omega )`$ is the response function of the P-I controller and $`\eta `$ is the response function of the VCO which is assumed to be frequency independent. The total response of the PLL depends crucially on the characteristic frequency of the P-I controller, $`f_{PI}`$, its gain $`P`$ and on the characteristic frequency of the tuning fork, $`f_G`$. Parameters are optimized, if

$$f_{PI}=f_G.$$
(16)
If $`f_{PI}<f_G`$ the response of the PLL is unnecessarily retarded while for $`f_{PI}>f_G`$ the response function develops an overshoot and tends to become unstable. If we regard $`V_{\mathrm{\Delta }f}`$ as the output and $`f_{res}`$ as the input of the PLL we find the solution of the linear response equations

$`V_f(\omega )=PLL(\omega )\left(f_{res}(\omega )f_0\right)`$

with the PLL-response function

$`PLL(\omega )={\displaystyle \frac{1}{\eta }}{\displaystyle \frac{L(\omega )}{1+L(\omega )}}`$

and the open loop response defined as

$`L(\omega ):=R(\omega )k(\omega )\eta .`$

The response becomes unstable at frequency $`\omega _0`$ if $`\text{Arg}[L(\omega _0)]=\pi `$ and $`\text{Abs}[L(\omega _0)]>1`$ (positive feedback). Stable feedback is therefore only achieved if $`P`$ is kept below a critical value $`P_c`$. In practice $`P_c`$ has to be determined experimentally.

At very low frequencies $`R(\omega )`$ is huge and $`L(\omega )1`$. In this limit $`PLL(\omega )`$ is constant with the value $`1/\eta `$. At frequencies much larger than $`f_{PI}=f_G`$ the function $`R(\omega )`$ is constant but $`k(\omega )`$ becomes very small such that $`L(\omega )1`$ and the response function decays like $`k(\omega )`$.

We now discuss the behaviour of this response function for the case that $`f_G=f_{PI}`$. In this case it reduces to

$`PLL(\omega )={\displaystyle \frac{1}{\eta }}{\displaystyle \frac{1}{1+i\frac{f}{Pa\eta f_G}}}={\displaystyle \frac{1}{\eta }}{\displaystyle \frac{1}{1+i\frac{f}{f_{PLL}}}},`$

where we have introduced $`a=X_{res}Q/(\pi f_0)`$, i.e. to a low-pass behaviour with the new characteristic frequency

$`f_{PLL}=Pa\eta f_G.`$

The product $`a\eta `$ acts like an additional $`P`$-factor and the effective $`P`$-factor is $`P_{eff}=Pa\eta `$. The bandwidth of the PLL is a factor of $`P_{eff}`$ larger than the bandwidth of the interaction sensor itself. This is shown in Fig. 4 as the curve labeled ‘PLL’ which was measured with $`\eta =0.1`$ Hz/V, $`a=71`$ V/Hz and $`P=2`$. Higher bandwidths can be achieved with higher $`P_{eff}`$, however, the stability condition $`P_{eff}<P_c`$ and noise considerations constitute the upper experimental limits.

### B Noise considerations

The dominating noise from the I-U converter can be represented as an equivalent noise source $`n_Y(\omega )`$ at the Y-output of the lockin. A linear response analysis starts from the equations

$`Y(\omega )`$ $`=`$ $`k(\omega )(f_{res}(\omega )f_d(\omega ))+n_Y(\omega )`$ (17)
$`V_f(\omega )`$ $`=`$ $`R(\omega )Y(\omega )`$ (18)
$`f_d(\omega )f_0`$ $`=`$ $`\eta V_f(\omega ).`$ (19)
The resulting noise term of the solution for $`f_G=f_{PI}`$ is

$$n_{PLL}(\omega )=P\frac{fif_G}{fif_{PLL}}n_Y(\omega ).$$
(20)
At frequencies $`ff_{PLL},f_G`$ the noise response is simply given by $`P`$, while for $`ff_{PLL},f_G`$ the noise is given by the ratio $`f_G/f_{PLL}=1/(Pa\eta )<1`$. The low frequency noise is reduced by the PLL. However, since $`f_{PLL}`$ is small compared to the lockin bandwidth of about 530 Hz, the low frequency noise suppression is not relevant in practice for the integrated noise spectrum and to a good approximation the output noise of the PLL $`\delta V_{\mathrm{\Delta }f}`$ is enhanced by the factor $`P`$ over the noise on the $`Y`$-output $`\delta Y`$:

$$\delta U_{\mathrm{\Delta }f}=P\delta Y.$$
(21)
A measurement of the PLL-noise is shown in Fig. 6. Here, $`\eta =0.1`$ Hz/V, the tuning fork excitation was 10 $`\mu `$V. The lockin time constant $`\tau =300`$ $`\mu `$s at 24 dB/oct is responsible for the strong cut-off of the noise spectrum above 500 Hz. The I-U converter output noise $`n_{IUC}`$ was amplified by a factor 1000 by the lock-in and therefore has a value of 800 $`\mu `$V/$`\sqrt{\text{Hz}}`$. In the frequency range between 5 and 500 Hz the PLL-noise is enhanced by the factor $`P=12`$. At the lowest frequencies the noise is reduced in agreement with the above discussion. The integrated noise spectrum corresponds to an effective frequency noise of about 20 mHz/$`\sqrt{\text{Hz}}`$. If desired, parameters can be chosen such that the frequency noise is suppressed well below 1 mHz/$`\sqrt{\text{Hz}}`$ at the cost of bandwidth. For scanning we typically work with a PLL-noise of 10 mHz/$`\sqrt{\text{Hz}}`$.

## VI The $`z`$-feedback

### A Linear response

The PLL-feedback is part of the $`z`$-feedback loop as shown in Fig. 7. The PLL-output voltage $`U_{\mathrm{\Delta }f}`$ is fed into the $`z`$-feedback P-I controller (PI) which drives the high-voltage amplifier (HV). This amplifier supplies the voltage for the $`z`$-piezo electrodes and thereby determines the tip-sample separation $`\mathrm{\Delta }z`$. As a result of the force-distance interaction characteristic, $`\mathrm{\Delta }z`$ and the tip-sample voltage $`U_{ts}`$ are translated into a resonance frequency of the interaction sensor \[see eq. (8)\] which is part of the PLL. For controlling the tip at constant height above the sample surface a certain $`U_{set}`$ is chosen corresponding to a certain frequency shift via $`f_{set}=\eta U_{set}`$. The feedback will then keep the tuning fork resonance frequency at $`f_{set}`$ by controlling $`\mathrm{\Delta }z`$ during a scan.

In the following we discuss the $`z`$-feedback in a similar fashion to the PLL. Here the relevant linear equations are

$`\mathrm{\Delta }f(\omega )`$ $`=`$ $`\alpha U_{ts}^{AC}(\omega )+\beta \mathrm{\Delta }z(\omega )`$ (22)
$`\mathrm{\Delta }z(\omega )`$ $`=`$ $`\mu \mathrm{\Delta }U_{HV}(\omega )`$ (23)
$`U_{\mathrm{\Delta }f}(\omega )`$ $`=`$ $`PLL(\omega )\mathrm{\Delta }f(\omega )`$ (24)
$`e(\omega )`$ $`=`$ $`U_{set}U_{\mathrm{\Delta }f}(\omega )`$ (25)
$`\mathrm{\Delta }U_{HV}(\omega )`$ $`=`$ $`PI(\omega )e(\omega )`$ (26)
Here $`\mu `$ is the (temperature-dependent) calibration of the $`z`$-piezo tube (e.g. 425 nm/70 V at 4.2 K), $`U_{set}`$ is the setpoint of the feedback, $`PI(\omega )`$ is the response of the $`z`$-feedback P-I controller and $`e(\omega )`$ is the error signal of the feedback loop. We are interested in the response of $`\mathrm{\Delta }z`$ as a result of surface roughness on the sample. We can simulate this roughness by applying a sinusoidal AC tip-sample voltage $`V_{ts}^{AC}`$. The solution is

$`\mathrm{\Delta }z=Z(\omega ){\displaystyle \frac{\alpha }{\beta }}U_{ts}^{AC}.`$

Here, the quantity $`\alpha /\beta U_{ts}^{AC}`$ plays the role of an effective surface roughness $`\mathrm{\Delta }z_{eff}`$ and

$`Z(\omega )={\displaystyle \frac{M(\omega )}{1+M(\omega )}}.`$

with the open loop response

$`M(\omega ):=\beta \mu PI(\omega )PLL(\omega ).`$

In the limit of low frequencies $`PLL(\omega )`$ is constant, $`PI(\omega )`$ dominates and $`\beta \mu PI(\omega )PLL(\omega )1`$ such that $`Z(\omega )`$ is constant and given by $`z_{eff}`$. In the limit of large frequencies this response decays like $`PLL(\omega )`$.

This feedback is stable, if $`\text{Abs}[M(\omega )]<1`$ at the frequency where $`\text{Arg}[M(\omega )]=\pi `$. This gives again an upper limit $`P_c^{(z)}`$ for possible $`P_z`$ values. $`P_c^{(z)}`$ depends on $`\beta `$ and therefore on the tip-sample interaction.

If we assume again that

$$f_{PLL}=f_z$$
(27)
the z-feedback response simplifies to

$`Z(\omega )={\displaystyle \frac{1}{1+if/f_{FB}}}`$

with

$$f_{FB}=\frac{P_z\beta \mu }{\eta }f_{PLL}=P_zP\beta \mu af_G.$$
(28)
Here the product $`\beta \mu a`$ acts as an additional $`P`$-factor and the effective value is $`P_{eff}^{(z)}=P_zP\beta \mu a`$ such that $`f_{FB}=P_{eff}^{(z)}f_G`$, i.e., the bandwidth of the $`z`$-feedback is again larger than the bandwidth of the PLL. This is demonstrated in Fig. 4 where the curve labeled ‘z-feedback’ was measured by varying the frequency of $`U_{ts}^{AC}(\omega )`$. The figure summarizes the resulting frequency characteristics of the two nested feedback-loops, i.e., the PLL and the $`z`$-feedback. It is demonstrated that the total system response is a low-pass, if all parameters are properly set. The bandwidth of the $`z`$-feedback can be increased by using higher $`P_{eff}`$, however, the stability criterion $`P_{eff}<P_c^{(z)}`$ and noise considerations constitute the experimental limits.

### B Noise considerations

If we consider the contribution of the dominant I-U-converter noise in the $`z`$-feedback the relevant equations are

$`\mathrm{\Delta }f(\omega )`$ $`=`$ $`\alpha U_{ts}^{AC}(\omega )+\beta \mathrm{\Delta }z(\omega )`$ (29)
$`\mathrm{\Delta }z(\omega )`$ $`=`$ $`\mu \mathrm{\Delta }U_{HV}(\omega )`$ (30)
$`U_{\mathrm{\Delta }f}(\omega )`$ $`=`$ $`PLL(\omega )\mathrm{\Delta }f(\omega )+n_{PLL}(\omega )`$ (31)
$`e(\omega )`$ $`=`$ $`U_{set}U_{\mathrm{\Delta }f}(\omega )`$ (32)
$`\mathrm{\Delta }U_{HV}(\omega )`$ $`=`$ $`PI(\omega )e(\omega )`$ (33)
The noise contribution $`n_{PLL}(\omega )`$ is the PLL output noise calculated earlier (see eq. (20)) which is dominated by the ouput noise of the I-U converter. The response to the PLL-output noise determined from these equations under the condition $`f_{PLL}=f_z`$ is

$`n_z(\omega )=\mu P_zP{\displaystyle \frac{fif_G}{fif_0}}n_Y(\omega ).`$

It is interesting to note that the bandwidth of the PLL does not enter the noise figure of the $`z`$-feedback. The integrated noise of the $`z`$-feedback is given by

$$\delta z=\mu P_zP\delta Y$$
(34)
Higher bandwidth of the $`z`$-feedback will lead to larger $`z`$-noise and for high quality imaging it is important to keep $`\delta z`$ as low as possible (typically of the order of 1 Å). On the other hand, reasonable scan speeds require high bandwidths of the order of several hundred Hz. This naturally leads to the question, how the optimum feedback parameters can be found from the above analysis for a given setup.

## VII Optimum feedback parameters

The optimum feedback parameters are given by the following conditions:

1. The characteristic frequency of the PLL-PI-controller is identical to the characteristic frequency $`f_G`$ of the (tuning fork) sensor.
2. The characteristic frequency of the $`z`$-feedback P-I controller is identical to the PLL-bandwidth.
3. The $`P_{eff}`$ parameter of the PLL is well below the critical value $`P_c`$ where the PLL-feedback becomes unstable.
4. The $`P_{eff}^{(z)}`$ parameter of the z-feedback is well below the critical value $`P_c^{(z)}`$ where the z-feedback becomes unstable.
5. The $`P_{eff}`$ parameter of the PLL is small enough to give acceptable frequency noise.
6. The $`P_{eff}^{(z)}`$ parameter of the $`z`$-feedback is small enough to give acceptable $`z`$-noise.

We repeat the corresponding equations (16), (27), (21) and (34) below:

$`f_{PI}`$ $`=`$ $`f_G={\displaystyle \frac{f_0}{2Q}}`$ (35)
$`f_z`$ $`=`$ $`f_{PLL}=P\eta {\displaystyle \frac{X_{res}}{2\pi }}`$ (36)
$`\delta U_{\mathrm{\Delta }f}`$ $`=`$ $`P\delta Y`$ (37)
$`\delta z`$ $`=`$ $`\mu P_zP\delta Y`$ (38)
From these equations the optimum feedback parameters can be found as follows: First, the tuning fork resonance curve is measured and $`f_G`$ is determined from its half width at the $`\sqrt{2}`$-maximum value. This frequency determines the characteristic frequency to be set on the P-I controller of the PLL via eq. (35). The gain $`P`$ of this P-I controller is set such that the output noise of the PLL \[$`\delta U_{\mathrm{\Delta }f}`$, eq. (37)\] corresponds to a reasonable frequency noise (typically 10 mHz for scanning in our setup). Next, the P-I parameters of the $`z`$-feedback are set such that eq. (36) is fulfilled. Finally, the gain $`P_z`$ of the $`z`$-feedback is adjusted according to the tolerable $`z`$-noise with the help of eq. (38). We typically keep the $`z`$-noise around 1 Å for normal scanning operation.

We emphasize that the optimum feedback parameters determined in this way do not depend on the tip shape or any details of the tip-sample interaction. Knowledge of the resonance curve of the sensor alone is sufficient for the optimum settings of the feedback.

For a given $`z`$-noise there is a certain freedom concerning the choice of the parameters $`P`$ and $`P_z`$, since according to eqs. (37) and (38) only their product needs to be kept constant. This freedom allows us to choose the bandwidth of the PLL either relatively high or relatively low in the span between $`f_G`$ and $`f_{FB}`$. If the PLL-bandwidth is chosen to be high, e.g. close to $`f_{FB}`$, the PLL will be able to track the resonance frequency with the same bandwidth as the $`z`$-feedback and $`U_{\mathrm{\Delta }f}`$ will be constant throughout a scan. However, if the PLL bandwidth is chosen to be low, e.g. closer to $`f_G`$, the PLL can not track the resonance frequency with the same bandwidth as the $`z`$-feedback and appreciable error signals arise whenever there is a step on the surface. This is demonstrated in Fig. 8. The two images in the top row are topographic images of a detail of a Hall bar structure with AFM-written oxide lines in the bottom right corner. The left (right) image was taken with small (large) PLL-bandwidth at fixed bandwidth of the $`z`$-feedback. The bottom row shows the corresponding frequency error signal. It can be seen that with a slow PLL, an additional image can be obtained in which the contours of edges appear very sharp. With a fast PLL this contrast disappears completely. This ‘differential’ image of the surface can be very useful since it eliminates sample tilt and amplifies edges of small height nicely.

## VIII Effect of Q and k on the bandwidth

It becomes clear from the above discussion that the bandwidth of the $`z`$-feedback is given by eq. (28)

$`f_{FB}=P_zP\beta \mu af_G=P_zP\beta \mu {\displaystyle \frac{X_{res}}{2\pi }}`$

i.e., it is independent of the $`Q`$-factor of the tuning fork cantilever. This is not obvious, since high $`Q`$ cantilevers have a low $`f_G`$ and are therefore known to slow down the response. However, implemented in our nested feedback system the quality factor of the sensor is canceled from the final expression for the bandwidth. The only impact on the performance stems from the fact that the sensor response can only be described as a low-pass, if changes in frequency $`\mathrm{\Delta }f`$ typically occurring during a scan are small compared to $`f_G=f_0/(2Q)`$.

The situation is different with the stiffness $`k`$ of the sensor which directly enters the quantity $`\beta `$, i.e., the slope of the $`\mathrm{\Delta }f(z)`$ curve at the operating point. Softer sensors will lead to larger values of $`\beta `$ and thereby increase the bandwidth of the $`z`$-feedback.

Similar arguments play a role for the comparison of scanning in the attractive or repulsive part of the $`\mathrm{\Delta }f(z)`$ curve. In the attractive part typical values of $`\beta `$ are much smaller than in the repulsive part. Therefore the bandwidth tends to be considerably larger for controlling in the repulsive regime.

## IX Comparison to phase control

We briefly want to discuss why the use of the PLL is advantageous compared to phase control, where the phase change measured as the $`Y`$-output of the lock-in in our case is directly used for the $`z`$-controller without the intermediate PLL.

From a similar linear response analysis as performed in this paper for our nested system one finds that the same bandwidths of the $`z`$-feedback can be achieved with pure phase control. Enhancement of the bandwidth for scanning operation is therefore not the striking argument for the use of the additional PLL feedback.

However, there are several reasons, why we use the PLL:

1. The PLL increases the bandwidth of the tuning fork response when the $`z`$-feedback is not yet controlling, e.g. during the approach of the tip. The PLL allows the use of reasonable approach speeds.
2. The PLL increases the dynamic range of the frequency detection. The Y-signal of the lock-in depends linearly on the frequency shift only in a limited range of frequencies around resonance. With high $`Q`$ sensors this range can become so narrow (100 mHz) that during tip approach the frequency shift runs out of this range. The PLL avoids this problem by tracking the resonance. The frequency range that can be used with PLL is mainly determined by the VCO-settings. We typically use a range of 2 Hz.
3. The PLL allows to image with an additional adjustable frequency contrast (see Fig. 8).
4. When one measures $`\mathrm{\Delta }f(z)`$-curves the PLL allows the clear distinction between conservative frequency shifts and dissipative effects .

## X Conclusion

In the first part of this paper we have shown that the tuning fork sensors that we use in our cryo-SFM setups can be well modeled by a single-harmonic oscillator despite the symmetry breaking wire that is attached to one tuning fork prong. In the second part of the paper we have analyzed the linear response of our nested feedback system comprising the interaction sensor, a phase-locked loop and a conventional $`z`$-feedback. We have demonstrated that the optimum feedback parameters of such a setup can be easily found and that they are independent of the tip shape and of the details of the tip-sample interaction. We have compared our setup to the conventional phase-control mode and discussed the advantages of the use of the phase-locked loop for high-Q cantilevers.


format: text

text:
# 1 Introduction

## 1 Introduction

Illuminating normally one of the two relevant faces of a wedge by a 10 GHz beam, Shelby et al.  observed that the angle of emergence of the beam from the other face changed sign when the teflon wedge was replaced by one made of the so–called left–handed material (LHM) . Assuming that their LHM is homogeneous and isotropic, and virtually disregarding dissipation, they concluded that it possesses a negative index of refraction.

Never mind that the so–called LHM does not actually possess handedness , and therefore is mis–named. Never mind that this material is actually a particulate composite  containing particles that are not electrically small — especially when the fiber–glass content of the LHM is factored in — and therefore is not homogenizable into a medium with local response properties . Never mind also that the LHM is actually anisotropic and must have some scattering losses.

A negative index of refraction is not needed per se, because all that it connotes is that the phase velocity is directed opposite to the velocity of energy transport  — as shown in the Appendix. But the phase velocity is an inessential concept , whether in vacuum or in materials. What really matter are the solutions of Maxwell equations that propagate energy on a desired trajectory, which is chosen epistemically.

Anyhow, virtually all publications on LHMs that I have come across ,  deal substantively only with non–dissipative materials, with just one exception . Perhaps, in time to come, samples of these materials with proven negligible dissipation at some frequency will emerge. In the meanwhile, they permit me to postulate an electromagnetic trinity comprising (i) vacuum, (ii) anti–vacuum, and (iii) nihility. This trinity illuminates the structure of electromagnetic theory, at least insofar as the relationship of phase velocity with Poynting vector is concerned.

## 2 Vacuum

Vacuum, of course, is well–known. Also called free space, it is a matter–free medium — the substrate on which microscopic electromagnetics is formulated and later homogenized for macroscopic research . The frequency–domain constitutive relations of vacuum are specified as

$$\begin{array}{c}𝐃(x,y,z,\omega )=ϵ_0𝐄(x,y,z,\omega )\hfill \\ 𝐁(x,y,z,\omega )=\mu _0𝐇(x,y,z,\omega )\hfill \end{array}\},$$
(1)
where $`x`$, $`y`$ and $`z`$ are the cartesian components of the position vector; $`\omega `$ is the angular frequency; $`ϵ_0=8.854\times 10^{12}`$ F m-1 and $`\mu _0=4\pi \times 10^7`$ H m-1. The phase velocity and the wavenumber of a plane wave in vacuum are co–parallel.

## 3 Anti–vacuum

The constitutive relations of anti–vacuum are postulated as

$$\begin{array}{c}𝐃(x,y,z,\omega )=ϵ_0𝐄(x,y,z,\omega )\hfill \\ 𝐁(x,y,z,\omega )=\mu _0𝐇(x,y,z,\omega )\hfill \end{array}\},$$
(2)
Among the more significant properties of anti–vacuum is that the phase velocity of a plane wave therein is opposite in direction to the wavevector. Hence, anti–vacuum is complementary to vacuum.

Clearly, anti–vacuum does not exist, but it could be simulated at a particular spatial frequency as follows: Suppose one could fabricate a LHM whose relative permittivity scalar $`ϵ_r(\omega )`$ is purely real–valued at a certain angular frequency $`\stackrel{~}{\omega }`$, and whose relative permittivity scalar $`\mu _r(\omega )`$ is such that $`\mu _r(\stackrel{~}{\omega })=ϵ_r(\stackrel{~}{\omega })`$. Let electrically small spheres of this LHM be homogeneously but randomly dispersed in vacuum. As both mediums would be impedance–matched, the effective constitutive parameters of the particulate composite thus formed would have to satisfy the relation $`\mu _r^{eff}(\stackrel{~}{\omega })=ϵ_r^{eff}(\stackrel{~}{\omega })`$ . Then, $`ϵ_r^{eff}(\stackrel{~}{\omega })=1`$ according to the Bruggeman formalism , provided the volume fraction of the LHM is

$$f_{LHM}=\frac{2}{3}\frac{ϵ_r(\stackrel{~}{\omega })2}{ϵ_r(\stackrel{~}{\omega })1},0f_{LHM}1.$$
(3)
As an example, $`f_{LHM}=0.8`$ if $`ϵ_r(\stackrel{~}{\omega })=4`$. Hence, this composite would mimic anti–vacuum at $`\omega =\stackrel{~}{\omega }`$.

## 4 Nihility

Nihility is the electromagnetic nilpotent, with the following constitutive relations

$$\begin{array}{c}𝐃(x,y,z,\omega )=\mathrm{𝟎}\hfill \\ 𝐁(x,y,z,\omega )=\mathrm{𝟎}\hfill \end{array}\}.$$
(4)
Wave propagation cannot occur in nihility, because $`\times 𝐄(x,y,z,\omega )=\mathrm{𝟎}`$ and $`\times 𝐇(x,y,z,\omega )=\mathrm{𝟎}`$ in the absence of sources therein . The directionality of the phase velocity relative to the wavevector in nihility is thus a non–issue.

The Bruggeman formalism is unable to predict the realization of nihility as a particulate composite either by mixing anti–vacuum and vacuum, or by mixing an LHM with an isoimpedant commonplace material (such as teflon). This is because the polarizability and the magnetizability (per unit volume) of an isotropic dielectric–magnetic sphere embedded in nihility do not depend on the constitutive parameters of the medium that the sphere is made of.

The Maxwell Garnett formalism , however, predicts that a homogeneous dispersion of electrically small, anti–vacuum spheres in vacuum is effectively equivalent to nihility, provided the volume fraction of anti–vacuum is $`0.25`$. Conversely, if matter were scooped out of anti–vacuum in the form of electrically small spheres, the resulting Swiss–cheese composite would also mimic nihility, if the volume fraction of anti–vacuum were $`0.75`$. Thus, nihility too could be mimicked by particulate composites at $`\omega =\stackrel{~}{\omega }`$.

## 5 Concluding Remarks

Whereas vacuum and anti–vacuum are mutually complementary in regard to the directionality of the phase velocity with respect to the Poynting vector, that issue does not arise for nihility because waves cannot propagate inside it. Of the three, vacuum alone is matter–free. Both anti–vacuum and nihility would have to be simulated — at a specific angular frequency — as particulate composites which would ultimately require the fabrication of virtually non–dissipative LHMs that are impedance–matched to vacuum. When would this become possible at all, let alone routinely, is anybody’s guess.

## Appendix

Let us consider an isotropic, dielectric–magnetic, homogeneous medium with relative permittivity $`ϵ_r=a\left(1i\mathrm{tan}\delta _ϵ\right)`$ and relative permittivity $`\mu _r=b\left(1i\mathrm{tan}\delta _\mu \right)`$ at a certain angular frequency $`\omega `$. Let $`a`$, $`b`$, $`m=+\sqrt{ab}`$, $`\mathrm{tan}\delta _ϵ`$ and $`\mathrm{tan}\delta _\mu `$ all be positive real, so that $`ϵ_r`$ and $`\mu _r`$ are in accord with the Lorentz model , with $`\omega `$ significantly larger than the resonant angular frequencies of the medium. The $`\mathrm{exp}\left(i\omega t\right)`$ time–dependence is implicit. To ensure that dissipation is moderate, let $`\mathrm{tan}\delta _ϵ<0.1`$ and $`\mathrm{tan}\delta _\mu <0.1`$.

The electromagnetic phasors of a planewave traveling along the $`+z`$ axis in the chosen medium are given by

$$\begin{array}{cc}𝐄(z,\omega )=A𝐮_x\mathrm{exp}\left(ik_0nz\right)\hfill & \\ 𝐇(z,\omega )=\frac{n}{\mu _r}\frac{A}{\eta _0}𝐮_y\mathrm{exp}\left(ik_0nz\right)\hfill & \end{array}\},$$
(5)
where $`A`$ is a complex–valued amplitude, $`k_0`$ is the vacuum wavenumber, $`\eta _0`$ is the intrinsic impedance of vacuum, and $`n^2=ϵ_r\mu _r`$. Writing $`n=n^{}+in^{\prime \prime }`$, we see that the time–averaged Poynting vector is given as

$$𝐏(z,\omega )=𝐮_z\mathrm{Re}\left\{\frac{n}{\mu _r}\right\}\frac{\left|A\right|^2}{2\eta _0}\mathrm{exp}\left(2k_0n^{\prime \prime }z\right)=𝐮_zP_z(z,\omega ),$$
(6)
the velocity of energy transport is co–parallel with $`𝐏(z,\omega )`$, and the phase velocity $`v_{ph}=\omega /k_0n^{}`$. Correct to first order in both $`\mathrm{tan}\delta _ϵ`$ and $`\mathrm{tan}\delta _\mu `$, we get $`n=\pm m\left[1i\left(\mathrm{tan}\delta _ϵ+\mathrm{tan}\delta _\mu \right)/2\right]`$ and $`\mathrm{Re}\left\{\frac{n}{\mu _r}\right\}=m/b`$.

For energy to flow along the $`+z`$ axis, we must have $`P_z>0`$. Furthermore, the directions of energy attenuation and energy flow must coincide (as befits a passive medium), so that $`n^{\prime \prime }>0`$. Both of these criterions lead to $`n^{}<0`$, which means that the phase velocity is pointed exactly opposite to the common direction of energy attenuation and flow.


text:
# Equilibration times in numerical simulation of structural glasses: Comparing parallel tempering and conventional molecular dynamics

## I MODELS AND DETAILS OF SIMULATION

The system we investigated is the monoatomic Lennard-Jones (LJ) model modified by adding a many-body anti-crystalline potential designed to inhibit crystallization angelani . The $`ϵ`$ and $`\sigma `$ parameter of the LJ potential are chosen as unit of energy and length respectively. The LJ potential is truncated and shift at $`2.5`$. The potential energy, which includes the anti-crystalline potential is

$`E(\stackrel{}{q}_m)`$ $`=V_{LJ}(\stackrel{}{q}_m)+`$ (3)
$`{\displaystyle \frac{1}{2}}\alpha {\displaystyle \underset{\stackrel{}{k}}{}}\theta (S_m(\stackrel{}{k})S_0)[S_m(\stackrel{}{k})S_0]^2.`$
where $`V_{LJ}`$ is the LJ part of the potential and the sum is over all $`\stackrel{}{k}`$ such that $`k_{max}\mathrm{\Delta }k<\stackrel{}{k}<k_{max}+\mathrm{\Delta }k`$. The other terms in Eq. set in only when any of the wavevector around the structure factor peak increases beyond the value $`S_0`$ and act by damping the unwanted crystalline like density fluctuation. We chose a number density $`\rho =1`$, and $`S_0=10`$, $`k_{max}=7.12`$, $`\alpha =0.83`$ and $`\mathrm{\Delta }k=0.34`$ for the anti-crystalline parameters as proposed in Ref.angelani . The integration timestep is $`0.0025`$, in time units of $`\sqrt{m\sigma ^2/ϵ}`$. The dynamics for this model has been previously studied newangelani . It has been shown that a fast decrease of the structural times takes place below $`T=1`$. The $`T`$ dependence of $`\tau `$ follows a power law in $`TT_x`$, with $`T_x0.475`$. $`T_x`$ has been identified with the ideal MCT for this model, an hypothesis supported also by an analysis of the $`T`$ dependence of the diffusive directions saddles .

The PT algorithm is identical to that encoded in Ref.KobYam  and we refer to that paper for details on the technique. The algorithm we implement uses $`M=14`$ identical non-interacting replicas each composed of $`N=256`$ particles. The $`14`$ temperatures are chosen to span an range from $`T=1.05`$ down to $`0.485`$, in particular the temperatures we used are the following: $`0.485`$, $`0.518`$, $`0.534`$, $`0.562`$, $`0.597`$, $`0.646`$, $`0.694`$, $`0.745`$, $`0.80`$, $`0.85`$, $`0.90`$, $`0.95`$, $`1.0`$, $`1.05`$.

The Hamiltonian of one replica is that of Eq. . All replicas evolve according to the standard Nose’ constant temperature MD simulation. Every $`1000`$ steps an attempt to exchange the scaling parameter of all pair of replicas with adjacent temperatures (swap of the $`\mathrm{\Lambda }`$ values) is performed using the following criterion: an exchange is accepted in Metropolis fashion, i.e. the acceptance ratio is:

$$w_{m,n}=\{\begin{array}{cc}1,\hfill & \mathrm{\Delta }_{m,n}0\hfill \\ \mathrm{exp}(\mathrm{\Delta }_{m,n}),\hfill & \mathrm{\Delta }_{m,n}>0\hfill \end{array}$$

where $`\mathrm{\Delta }_{m,n}=\beta _0(\mathrm{\Lambda }_n\mathrm{\Lambda }_m)[E(q_m)E(q_n)]`$. The events with $`i=0,2,4,\mathrm{}`$ or $`i=1,3,5,\mathrm{}`$ are repeated alternatively every $`1000`$ integration steps.

The outcome of such calculation are, in principle, equilibrium configurations in the canonical ensemble at the $`M`$ different temperatures.

To estimate the time requested to the PT algorithm to equilibrated all replicas we start our PT algorithm with $`M`$ replicas extracted from a previously generated ensemble of equilibrium configurations at $`T=1.05`$. At this $`T`$, the structural relaxation time is of the order of 1000 steps and hence generation of equilibrium configurations with conventional MD does not pose any problem. By starting with this ensemble of configurations, the PT equilibration time is by construction $`0`$ for the highest temperature. We performed $`20`$ of such independent PT runs to improve the statistic. Each of such run lasts two million time step. Hence, in the PT part of the work, the equation of motion have been integrated more than 40 million times.

The same starting configurations ($`T=1.05`$) are also used as initial configurations for conventional constant temperature MD simulations at several bath temperatures $`T`$, to compare the rate of equilibration of PT and conventional MD algorithms. For each temperature of these MD simulations we performed $`16`$ independent runs to improve statistics.

Local minima configurations have been calculated via conjugate gradient minimization. The minimization process is considered completed when the potential energy change associated with one iteration is less than $`10^{15}`$ to ensure a great accuracy.

## II RESULTS

We focus here on the evolution in time of the inherent structure energy $`e_{IS}`$, comparing the PT and MD procedures. Recent work epl ; mossacondmat ; jpcm2001  has provided evidence that following a $`T`$ change, the system response is characterized by two different time scales. A short time scale, related to the equilibration of the system within a well defined basin of the energy landscape and a slow time scale related to the search for basin of the “right” depth. The evolution of one-time quantities carries on information on these two time scales. After a sudden change of $`T`$, a very fast decreases of the vibrational energy — corresponding to the fast equilibration to the new bath $`T`$ of the intra-basin vibrational motions — is observed. This fast change is followed by a much slower decrease which corresponds to the slow decrease of the basin’s depth. The absolute change in $`e_{IS}`$ during the aging process is significantly smaller than the change in the total potential energy and hence it requires a careful analysis to be detected. A way which has been proven fruitful to separate large fast component and small slow component is to monitor directly the evolution of $`e_{IS}`$. Building on the expertise developed in recent years, we adopt this indicator as effective tool to monitor the equilibration of the system in configuration space.

We note on passing that, since $`e_{IS}`$ is a small component of the potential energy — being the intra-basin vibrational part dominating — a nice scaling of the total potential energy distribution may not guarantee perfect equilibration.

Figure 1 compares the time evolution of the inherent structure energy in the conventional MD (Top) and in the PT (Bottom) runs. In both cases, by construction, the initial $`e_{IS}`$ coincides with the equilibrium value of $`e_{IS}`$ at $`T=1.05`$. As time goes on, each replica starts to explore larger and larger parts of the configuration space selecting configuration with lower and lower values of $`e_{IS}`$. The equilibration process lasts until the equilibrium value of $`e_{IS}`$ is reached. The same picture applies to the conventional MD case, where the simulation indeed reproduces the aging process following a $`T`$ jump from $`T=1.05`$ to the new bath temperature.

Fig. 2 compares the time evolution of $`e_{IS}`$ for PT and MD simulations at three different temperatures. In all cases, we find clear indications that the equilibration of the slow degrees of freedom does not depend on the procedure adopted.

## III DISCUSSION

The data shown in Fig.2 very clearly show that for the Lennard Jones case investigated in this article no improvement in equilibration rates is achieved by implementing the PT algorithm.

For the related model of binary mixture Lennard Jones kobrev , the number of basins with depth $`e_{IS}`$, in the range of $`e_{IS}`$ values characteristic of the PES region explored above $`T_x`$ is well represented by a gaussian distribution. The total number of basins has been shown to scale with the size $`N`$ of the system as $`e^{\alpha N}`$ with $`\alpha 0.8`$ sciortinoentropyprl . For our 256 atom case, this correspond to about $`10^{87}`$ basins. Under such complicated potential energy landscape conditions, an unbiased search for the location of the deepest basin would require an order of $`10^{87}`$ attempts! In this respect, it is possible that the rate of equilibration at low temperatures is significantly controlled by a simple entropic effect. This could explain why the possibility of overcoming barriers with higher probability offered by PT does not favor a faster equilibration process. This picture is also consistent with the fact that in the $`T`$-region explored (which is still above $`T_x`$) saddle dominated dynamics is dominant. Recent instantaneous normal modes analysis normalmodes  has indeed provided evidence that above $`T_x`$ the system explore mostly regions of the potential energy landscape which are characterized by a large number of negative curvature directions. No activated processes are requested in this condition to change local basin.

It would be interesting to find out if the PT algorithm may be valuable in studying strong liquids, for which less relevant changes in the PES are taking place on cooling as compared to fragile liquids saikevoivod-nature  and for which activated processes are dominant at low $`T`$. Preliminary indications kobprivatecom  seems to suggest that this may be the case. It would also be important to correlate the efficiency of PT with the structure of configuration space and connectivity between distinct potential energy surface basins. A possible line of research could be to compare $`e_{IS}(t)`$ for PT and MD in clusters with different disconnectivity graph disconnectivity  types.

The fact that in the explored $`T`$-region equilibration times are not improved in the PT case is consistent with recent finding in the field of protein simulations. This analogy is not at all surprising due to the significant similarities in the problem of glass formation and protein folding walesscience . Again, the similarity points to the structure of the PES as key element in the control of the characteristic times.

To conclude, we like to call the reader attention on the fact that to study a fixed $`T`$-range, the PT technique requires a number of replicas which increases linearly with the system size, to guarantee proper overlap in the potential energy distributions of adjacent replicas and hence a significant replica exchange rate. Moreover, all replicas have to be simulated for the same total time interval. This time is fixed by the lowest temperature, which is characterized by a relaxation time which may well be several order of magnitude smaller than the one of the top temperature. Both effects concur in making PT a not very convenient algorithm for simulating structural glasses as compared to conventional MD. Indeed, in MD, the total simulation time at each temperature can be chosen to scale with the structural relaxation time. Of course, bookkeeping facilities are enhanced in the PT case.

###### Acknowledgements.

The authors warmly thank Walter Kob for useful discussions. We acknowledge support from INFM - Iniziativa Calcolo Parallelo and INFM PRA HOP.


text:
# Noncommutativity and Model Building

## Acknowledgements

We would like to thank Giovanni Amelino-Camelia, Harald Dorn, Gordy Kane and Tim Morris for discussions. This work was partially supported by a PPARC SPG grant.