text:
# Thermal and Nonthermal Pion Enhancements with Chiral Symmetry Restoration

## 1 Introduction

It is generally believed that there are two QCD phase transitions in hot and dense nuclear matter. One of them is related to the deconfinement process in moving from a hadron gas to a quark-gluon plasma, and the other one describes the transition from the chiral symmetry breaking phase to the phase in which it is restored. Since the deconfinement and restored phases are only intermediate states in relativistic heavy ion collisions which are considered to be the way to generate the two transitions in laboratories, the key problem is how to extract possible signatures of the phase transitions hidden in the final state distributions.

One way to study the signatures of the chiral transition is to compare the properties of hadronic systems with and without undergoing a chirally symmetric phase. In recent years many such signatures have been proposed, such as the excess of pions due to a rapid thermal and chemical equilibration in the symmetric phase, the excess of low energy photon pairs by pion annihilation in hot and dense medium, the dilepton enhancement from leptonic decay of sigma, and the enhancement of the continuum threshold in the scalar channel. Most of the signatures are associated with the changes of sigma properties at finite temperatures and densities. In the vacuum there is no obvious reason for the presence of the scalar meson sigma in the study of chiral properties. However, there is a definite need for sigma at finite temperatures and densities in the investigation of chiral symmetry restoration. With increasing $`T`$ and $`\mu `$ sigma changes its character from a resonance with large mass to a bound state with small mass. This medium dependence of sigma mass makes important sense not only in the above mentioned particle yields, but also in thermodynamics. The contribution from sigma to any thermodynamic function can for most case be neglected, because it is heavy, while it plays an important role in the region around the critical point of the chiral phase transition, since there $`m_\sigma =0`$.

From the recent QCD model calculations, the chiral phase transition in the chiral limit with massless pions is of second order at high temperatures, and of first order at high densities. Therefore, there is a tricritical point $`P`$ in the $`T\mu `$ plane which separates the first- and second-order phase transition. In the real world with nonzero pion mass, the second order phase transition becomes a smooth crossover and the tricritical point $`P`$ becomes an endpoint $`E`$ of a first-order phase transition line. In relativistic heavy ion collisions, for the choice of control parameters like colliding energy and impact parameter such that the freeze-out of the system occurs near the point $`E`$, sigma is one of the most numerous species at the freeze-out, since it is nearly massless. These light sigmas at the freeze-out will lead to large event-by-event fluctuations of final state pions due to the $`\sigma `$-exchange in the pion thermodynamic potential. Another direct signature of this critical endpoint $`E`$ is a nonthermal excess of pions from the hadronic decay of sigma. With the expansion of the hadronic system produced in relativistic heavy ion collisions, the in-medium sigma mass rises towards its vacuum value and eventually exceeds the $`\pi \pi `$ threshold. As the $`\sigma \pi \pi `$ coupling is large, the decay proceeds rapidly. Since this process occurs after freeze-out, the pions generated by it do not have a chance to thermalize. Thus, it is expected that the resulting pion spectrum will have a nonthermal enhancement at low transverse momentum.

The present paper is devoted to a systematic and detailed investigation of the above idea on excess pions produced by sigma decay in a chiral quark model in the whole $`T\mu `$ plane, without the assumption that the freeze-out occurs at the chiral transition. In the chiral limit, we know very well that the sigma decay into two pions starts at the phase transition point, and the freeze-out happens later. Therefore, the sigma decay before and after the freeze-out will result in thermal and nonthermal enhancements of pions, respectively. Our motivation here is to see that in which temperature and density region of freeze-out there is a strong nonthermal enhancement.

We proceed as follows. In Section 2 we calculate the sigma decay rate as a function of $`T`$ and $`\mu `$ to the lowest order in a $`1/N_c`$ expansion in the framework of the Nambu-Jona-Lasinio model (NJL) at quark level, and discuss the importance of the Bose-Einstein statistics at high $`T`$ and the dominance of the first-order transition at high $`\mu `$. The in-medium pion mass drops down with decreasing $`T`$ and $`\mu `$. When it approaches to its vacuum value, the freeze-out happens. With this definition of freeze-out and the help of scaling hydrodynamics to describe the expansion of the hadronic system, we study in Section 3 the thermal and nonthermal pion enhancements due to $`\sigma `$ decay. Finally we give our summary and point out open problems that remain.

## 2  sigma decay rate at finite temperatures and densities

One of the models that enables us to see directly how the dynamical mechanisms of chiral symmetry breaking and restoration operate is the NJL model applied to quarks. Within this model one can obtain the hadronic mass spectrum and the static properties of mesons remarkably well. In particular, one can recover the Goldstone mode, and some important low-energy properties of current algebra such as the Goldberger-Treiman and GellMann-Oakes-Renner relations. Because of the contact interactions between quarks that is introduced in this model, there is of course no confinement. A further consequence of this feature is that the model is non-renormalizable, and it is necessary to introduce a regulator $`\mathrm{\Lambda }`$ that serves as a length scale in the problem. The sigma decay rate at zero density in this model was firstly studied by Hatsuda and Kunihiro. In order to investigate the thermal and nonthermal pion yields by sigma decay, we recalculate here the decay rate in the whole $`T\mu `$ plane, and analyze its relation with the first-order transition in high density region.

The two-flavor version of the NJL model is defined through the lagrangian density,

$$L_{NJL}=\overline{\psi }(i\gamma ^\mu _\mu m_0)\psi +G[(\overline{\psi }\psi )^2+(\overline{\psi }i\gamma _5\tau \psi )^2],$$
(1)
where only the scalar and pseudoscalar interactions corresponding to $`\sigma `$ and $`\pi `$ mesons, respectively, are considered, $`\psi `$ and $`\overline{\psi }`$ are the quark fields, the operate $`\tau `$ is the $`SU(2)`$ isospin generator, $`G`$ is the coupling constant with dimension $`GeV^2`$, and $`m_0`$ is the current quark mass.

In an expansion in the inverse number of colors, $`1/N_c`$, the zeroth order (Hartree) approximation for quarks together with the first-order (RPA) approximation for mesons gives a self-consistent treatment of the quark-meson plasma in the NJL model. In the chiral limit the tricritical point $`P`$ is located at $`T_P=79MeV,\mu _P=280MeV`$ in the $`T\mu `$ plane. In the case of explicit chiral symmetry breaking the endpoint of the first-order transition line is shifted with respect to the tricritical point $`P`$ towards large $`\mu `$, its position is $`T_E=23MeV,\mu _E=325MeV`$. To the lowest order, namely $`\left(1/N_c\right)^{1/2}`$, the Feynman diagram for the decay process $`\sigma 2\pi `$ is sketched in Fig.1. The partial decay rate in the rest frame of sigma is given in terms of the Lorentz-invariant matrix element $`M`$ by

$$\frac{d\mathrm{\Gamma }_{\sigma 2\pi }}{d\mathrm{\Omega }}=\frac{1}{32\pi ^2}\frac{|𝐩|}{m_\sigma ^2}|M|^2,$$
(2)
where $`|𝐩|=\sqrt{\frac{m_\sigma ^2}{4}m_\pi ^2}`$ is the pion momentum, $`m_\sigma `$ and $`m_\pi `$ are respectively $`\sigma `$\- and $`\pi `$-mass. Using an obvious notation for the $`M`$-matrix element, one has

$$M=g_\sigma g_\pi ^2A_{\sigma \pi \pi }$$
(3)
with

$$A_{\sigma \pi \pi }=T\underset{n}{}e^{i\omega _n\eta }\frac{d^3𝐪}{(2\pi )^3}TrS(𝐪,i\omega _n)\mathrm{\Gamma }_\pi S(𝐪+𝐩,i\omega _n+\frac{m_\sigma }{2})\mathrm{\Gamma }_\pi S(𝐪,i\omega _n+m_\sigma ).$$
(4)
Here $`\omega _n=(2n+1)\pi T,n=0,\pm 1,\pm 2,\mathrm{}`$ are fermionic Matsubara frequencies, $`\mathrm{\Gamma }_\pi =i\gamma _5\tau `$ is the pseudoscalar interaction vertex. The quark propagator is denoted by $`S(𝐪,i\omega _n)=\left(\gamma _0(i\omega _n+\mu )\gamma 𝐪+m_q\right)/\left((i\omega _n+\mu )^2E_q^2\right)`$ with the quark energy $`E_q=\sqrt{m_q^2+𝐪^2}`$. In Eq.(4), $`Tr`$ refers to the trace over color, flavor and spinor indices. The pion-quark-quark and sigma-quark-quark coupling strengths in Eq.(3) are determined in the model via

$`g_\pi ^2(T,\mu )={\displaystyle \frac{\mathrm{\Pi }_\pi (k_0,\mathrm{𝟎};T,\mu )}{k_0^2}}|_{k_0^2=m_\pi ^2},`$
$`g_\sigma ^2(T,\mu )={\displaystyle \frac{\mathrm{\Pi }_\sigma (k_0,\mathrm{𝟎};T,\mu )}{k_0^2}}|_{k_0^2=m_\sigma ^2},`$ (5)
where $`\mathrm{\Pi }_\pi `$ and $`\mathrm{\Pi }_\sigma `$ are the standard mesonic polarization functions for pion and sigma. The dynamically generated quark mass $`m_q(T,\mu )`$ and the meson masses $`m_\pi (T,\mu )`$ and $`m_\sigma (T,\mu )`$ are calculated using the usual gap equations in the Hartree approximation for quarks and RPA approximation for mesons. After evaluation of the Matsubara sum in Eq.(4), one has

$`A_{\sigma \pi \pi }(T,\mu )=`$ $`4m_qN_cN_f{\displaystyle }{\displaystyle \frac{d^3𝐪}{(2\pi )^3}}{\displaystyle \frac{f_F(E_q\mu )f_F(E_q+\mu )}{2E_q}}\times `$ (6)
$`{\displaystyle \frac{8(𝐪𝐩)^2(2m_\sigma ^2+4m_\pi ^2)𝐪𝐩+m_\sigma ^4/22m_\sigma ^2E_q^2}{(m_\sigma ^24E_q^2)\left((m_\pi ^22𝐪𝐩)^2m_\sigma ^2E_q^2\right)}},`$
where $`f_F(x)=\left(1+e^{x/T}\right)^1`$ is the Fermi-Dirac distribution function.

Sigma can decay into neutral and charged pions. By considering the exchange contribution for $`\sigma 2\pi _0`$ due to the identical particle effect, the total decay rate can be written as

$`\mathrm{\Gamma }_{\sigma 2\pi }(T,\mu )`$ $`=\mathrm{\Gamma }_{\sigma 2\pi _0}(T,\mu )+\mathrm{\Gamma }_{\sigma \pi _+\pi _{}}(T,\mu )`$ (7)
$`={\displaystyle \frac{3}{8\pi }}{\displaystyle \frac{\sqrt{\frac{m_\sigma ^2}{4}m_\pi ^2}}{m_\sigma ^2}}g_\sigma ^2g_\pi ^4|A_{\sigma \pi \pi }(T,\mu )|^2\left(1+2f_B({\displaystyle \frac{m_\sigma }{2}})\right),`$
where we have taken into account the Bose-Einstein statistics in terms of the distribution function $`f_B(x)=\left(e^{x/T}1\right)^1`$ for the final state pions.

From the comparison of $`\sigma \pi \pi `$ coupling strength $`g_{\sigma \pi \pi }`$ defined through the lagrangian $`L_Ig_{\sigma \pi \pi }\sigma \pi \pi `$ and the $`M`$-matrix element (3) in the NJL model, we have

$$g_{\sigma \pi \pi }(T,\mu )=2M(T,\mu )=2g_\sigma g_\pi ^2A(T,\mu ).$$
(8)
In the numerical calculations we used the dynamical quark mass $`m_q=0.32GeV`$, the pion decay constant $`f_\pi =0.093GeV`$ and the pion mass $`m_\pi =0.134GeV`$ in the vacuum ($`T=\mu =0`$) as input to determine the $`3`$ parameters in the NJL model, namely the coupling constant $`G=4.93GeV^2`$, the momentum cutoff $`\mathrm{\Lambda }=0.653GeV`$ and the current quark mass $`m_0=0.005GeV`$.

We first discuss the temperature dependence of $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$ in the case of zero baryon density which corresponds to the central region of relativistic heavy ion collisions. The coupling strength $`g_{\sigma \pi \pi }`$ is about $`2GeV`$ in the vacuum and almost a constant in the temperature region $`T<0.1GeV`$, and then decreases smoothly with increasing temperature, as shown in Fig.2. The energy condition for the decay process $`\sigma 2\pi `$ to happen is $`m_\sigma 2m_\pi `$ in the rest frame of $`\sigma `$. At the threshold temperature $`T0.19GeV`$ determined by $`m_\sigma (T)=2m_\pi (T)`$, $`g_{\sigma \pi \pi }`$ has a sudden jump from about $`0.3/fm`$ to zero. In a wide temperature region, the decay rate (7) goes up gradually from its vacuum value $`0.45/fm`$ to the maximum $`0.5/fm`$ at $`T0.17GeV`$, then it drops down rapidly but continuously and becomes zero at the threshold temperature. The fact that the maximum decay rate is not in the vacuum where the mass difference between $`\sigma `$ and $`\pi `$ has the largest value, but at $`T0.17GeV`$ is beyond our expectation. It originates from the important Bose-Einstein statistics of the final state pions at high temperatures. This can be seen clearly in Fig.2 where we show the decay rates with and without consideration of the Bose-Einstein statistics. The rapid decrease of $`\mathrm{\Gamma }_{\sigma 2\pi }`$ at high temperatures is due to the behavior of $`g_{\sigma \pi \pi }`$ and the phase space factor $`\sqrt{m_\sigma ^2/4m_\pi ^2}/m_\sigma ^2`$.

The density dependence of $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$ is shown in Fig.3. Since the threshold temperature for $`\sigma `$ decay decreases with increasing chemical potential, the crucial Bose-Einstein statistics effect at high temperatures is gradually washed out as the chemical potential goes up. There is almost no difference between with and without such statistics when $`\mu 0.25`$. Instead of the statistics the first-order chiral transition dominates the behavior of $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$ at high densities. When the system passes through the first-order transition line from the chirally symmetric phase to the broken phase, all the masses, $`m_q,m_\pi `$ and $`m_\sigma `$ are discontinuous, $`m_q`$ and $`m_\sigma `$ jump up and $`m_\pi `$ jumps down. Therefore, the mass difference between $`\sigma `$ and $`\pi `$ has a jump from less than $`2m_\pi `$ to larger than $`2m_\pi `$, and the $`\sigma `$ decay starts not at the threshold, but at some value beyond the threshold. This means that the decay in the first-order transition region happens suddenly with a nonzero rate, the temperature dependence of the decay rate behaves like a step function. This is shown in Fig.3 for $`\mu =0.33GeV\mu _E`$. When approaching to the endpoint $`E`$ from the first-order transition side, the jump of the mass difference disappear. The decay rate becomes continuous, while the coupling strength has still a jump.

To see the parameter dependence of the above discussion, we considered another set of parameters, $`G=5.46GeV^2,\mathrm{\Lambda }=0.632GeV`$ and $`m_0=0.0055GeV`$. It makes shifts of the magnitudes of $`g_{\sigma \pi \pi }`$ and $`\mathrm{\Gamma }_{\sigma 2\pi }`$, but does not change their shapes, especially the behavior around the threshold temperature or the chiral phase transition point.

Let’s summarize the calculations of the $`\sigma `$ decay rate. When the chemical potential is extremely high, the first-order chiral transition results in an almost temperature-independent decay rate. For most of the decay processes in the crossover region, the maximum decay rate does not locate in the vacuum, but is very close to the threshold point. Only when the chemical potential is around $`0.25GeV`$ which is close to the endpoint $`E`$, the decay rate increases continuously and not very fast in the beginning, and then it behaves like a constant. As we discussed in the introduction, the freeze-out of the pions may occur after $`\sigma `$ decay. According to the decay rate shown in Fig.3, when the decay happens in the first-order transition region or most of the places in the crossover region in the $`T\mu `$ plane, a large fraction of $`\sigma `$’s may decay into pions before the freeze-out, and these pions are thus thermalized before their freeze-out, only those pions generated after the freeze-out contribute to the nonthermal enhancement. However, when the $`\sigma `$ decay happens with $`\mu 0.25GeV`$, most of the pions may be produced after the freeze-out, and therefore, there may be a strong nonthermal enhancement in this case. In the next section we discuss the nonthermal enhancement quantitatively.

## 3  thermal and nonthermal pion enhancements

Before we calculate the pion numbers with the $`\sigma `$ decay rate as a function of temperature $`T`$ and chemical potential $`\mu `$ given in the last section, we need to treat $`T,\mu `$ as time-dependent. In the baryon-free ($`\mu =0`$) region of relativistic heavy ion collisions, the Bjorken’s scaling hydrodynamics is often used to describe the space-time evolution of hadronic thermodynamics. The time dependence of temperature in this scenario is given by

$$T=T_d\left(\frac{t_d}{t}\right)^\alpha $$
(9)
with $`\alpha =1/3`$. Here we have defined $`t_d`$ to be the time at which $`\sigma `$ decay begins and the threshold temperature $`T_d`$ satisfies the mass relation $`m_\sigma (T_d)=2m_\pi (T_d)`$. We will be interested in times $`tt_d`$. In hydrodynamic models the initial time corresponding to the highest temperature of the fluid is normally taken to be $`1fm`$. Since the sigma decay starts in the early stage of the fluid, we choose $`t_d=2fm`$ in our numerical calculations.

For heavy ion collisions with full nuclear stopping, the lifetime of the baryon-rich fluid is much shorter compared with that of the baryon-free fluid, the space evolution can be considered as a homogeneous 3-dimensional expansion (9) with $`\alpha =1`$ instead of $`\alpha =1/3`$. The real relativistic heavy ion collisions are between these two limits. We still use the expansion mechanism(9) as the time dependence of temperature for any chemical potential $`\mu `$, but the exponent $`\alpha `$ is now $`\mu `$-dependent and parameterized by

$$\alpha (\mu )=\left(1+2\frac{\mu }{\mu _{max}}\right)\alpha (0)$$
(10)
with $`\alpha (0)=1/3`$. Here the maximum chemical potential $`\mu _{max}`$ corresponding to the full stopping limit is chosen as $`0.4GeV`$ in the NJL model.

With the known decay rate $`\mathrm{\Gamma }(T,\mu )`$ and the relation between temperature and time for any chemical potential, the number of $`\sigma `$’s present at time $`t`$ is determined by

$$\frac{dN_\sigma (t)}{dt}=\mathrm{\Gamma }(t)N_\sigma (t),$$
(11)
and is therefore

$$N_\sigma (t)=N_\sigma (t_d)e^{_{t_d}^t\mathrm{\Gamma }(t^{})𝑑t^{}},$$
(12)
where $`N_\sigma (t_d)`$ is the maximum number at the beginning time $`t_d`$,

$`N_\sigma (t_d)=V(t_d)n_\sigma (t_d),`$
$`n_\sigma (t_d)={\displaystyle \frac{d^3𝐩}{(2\pi )^3}\frac{1}{e^{ϵ_\sigma (T_d)/T_d}1}}`$ (13)
with the $`\sigma `$ energy $`ϵ_\sigma =\sqrt{m_\sigma ^2+𝐩^2}`$ and the volume $`V(t_d)`$ of the system at time $`t_d`$.

Since each sigma yields two pions, the number of pions generated by $`\sigma `$ decay at time $`t`$ is related to the number of sigma by

$$N_\pi (t)=2\left(N_\sigma (t_d)N_\sigma (t)\right).$$
(14)
As the time $`t\mathrm{}`$, $`N_\sigma (t\mathrm{})=0`$, sigma disappears, and the pion number reaches its maximum $`N_\pi (t\mathrm{})=2N_\sigma (t_d)`$.

The conventional definition for freeze-out is that the time between two successful scattering events is larger than the collective expansion time scale,

$$\tau _{scatt}\tau _{exp}.$$
(15)
$`\tau _{scatt}`$ can be calculated from the thermally averaged $`\pi \pi `$ cross section in the NJL model, and $`\tau _{exp}=1/_\mu u^\mu `$ is simply reduced to $`\tau _{exp}=t`$ in the scaling hydrodynamics. Thus the freeze-out time $`t_f`$ and the freeze-out temperature $`T_f`$ are determined by $`\tau _{scatt}=t_f`$ and the expansion mechanism (9). With the two temperature scales $`T_d`$ and $`T_f`$ we can divide the total pions generated by $`\sigma `$ decay into thermal and nonthermal pions. For $`T_d>T_f`$, the produced pions in the interval $`T_f<T<T_d`$ have time to thermalize before the freeze-out, they lead to an enhancement of thermal pions, while in the interval $`T<T_f`$ the produced pions do not get a chance to thermalize, they result in a nonthermal enhancement of pions with low momentum. For $`T_d<T_f`$, $`\sigma `$ decay yields no thermal pions, all the resulted pions are nonthermal pions.

The thermal and nonthermal enhancements are defined as ratios of the numbers of the pions produced by $`\sigma `$ decay before and after the freeze-out to the number of the pions produced directly by the freeze-out of the original thermal pion gas,

$`r_{th}(t)={\displaystyle \frac{N_\pi (t)}{N_\pi ^{dir}(t_f)}},t_dtt_f,`$
$`r_{nt}(t)={\displaystyle \frac{N_\pi (t)N_\pi (t_f)}{N_\pi ^{dir}(t_f)}},t_ft`$ (16)
where $`N_\pi ^{dir}(t_f)`$ is determined by the freeze-out temperature $`T_f`$, the volume of the system at time $`t_f`$, and the pion energy $`ϵ_\pi =\sqrt{m_\pi ^2+𝐩^2}`$,

$`N_\pi ^{dir}(t_f)=V(t_f)n_\pi ^{dir}(t_f),`$
$`n_\pi ^{dir}(t_f)={\displaystyle \frac{d^3𝐩}{(2\pi )^3}\frac{3}{e^{ϵ_\pi (T_f)/T_f}1}}.`$ (17)
Obviously the measurable thermal and nonthermal enhancements are $`r_{th}(t_f)`$ and $`r_{nt}(t\mathrm{})`$.

With the in-medium pion momentum $`p=\sqrt{\frac{m_\sigma ^2}{4}m_\pi ^2}`$, the expansion mechanism (9), and the pion distribution (14), we determine the normalized momentum distribution of the nonthermal pions,

$$\frac{1}{N_{nt}}\frac{dN_{nt}}{dp}=\frac{e^{_{t_d}^t\mathrm{\Gamma }(t^{})𝑑t^{}}\mathrm{\Gamma }(t)\theta (tt_f)\frac{dt}{dp}}{_{t_f}^{\mathrm{}}e^{_{t_d}^t\mathrm{\Gamma }(t^{})𝑑t^{}}\mathrm{\Gamma }(t)𝑑t},$$
(18)
and the observable mean momentum of the nonthermal pions,

$$\overline{p}_{nt}=\frac{p\frac{dN_{nt}}{dp}𝑑p}{\frac{dN_{nt}}{dp}𝑑p}=\frac{_{t_f}^{\mathrm{}}p(t)e^{_{t_d}^t\mathrm{\Gamma }(t^{})𝑑t^{}}\mathrm{\Gamma }(t)𝑑t}{_{t_f}^{\mathrm{}}e^{_{t_d}^t\mathrm{\Gamma }(t^{})𝑑t^{}}\mathrm{\Gamma }(t)𝑑t}.$$
(19)
Since it is difficult to separate out the nonthermal pions experimentally, we calculate also the normalized momentum spectra of the total pions including the direct pions, the thermal pions from sigma decay, and the nonthermal pions from sigma decay,

$`{\displaystyle \frac{1}{N_{total}}}{\displaystyle \frac{dN_{total}}{dp}}={\displaystyle \frac{1}{n_\pi ^{dir}(t_f)+2\frac{t_d}{t_f}n_\sigma (t_d)}}\times `$ (20)
$`\left(\left(1+2{\displaystyle \frac{t_d}{t_f}}{\displaystyle \frac{n_\sigma (t_d)}{n_\pi ^{dir}(t_f)}}\left(1e^{_{t_d}^{t_f}\mathrm{\Gamma }(t^{})𝑑t^{}}\right)\right){\displaystyle \frac{dn_\pi ^{dir}(t_f,p)}{dp}}+2{\displaystyle \frac{t_d}{t_f}}n_\sigma (t_d)e^{_{t_d}^t\mathrm{\Gamma }(t^{})𝑑t^{}}\mathrm{\Gamma }(t){\displaystyle \frac{dt}{dp}}\theta (tt_f)\right),`$
where we have used

$$\frac{V(t_d)}{V(t_f)}=\frac{t_d}{t_f}$$
(21)
in the scaling hydrodynamics, and $`\frac{dn_\pi ^{dir}(t_f,p)}{dp}`$ is defined as

$$\frac{dn_\pi ^{dir}(t_f,p)}{dp}=\frac{p^2}{2\pi ^2}\frac{3}{e^{ϵ_\pi (T_f)/T_f}1}.$$
(22)
We now discuss our numerical calculations about the thermal and nonthermal pion distributions. Fig.4 shows the sigma number and pion number scaled by the corresponding maximum value,

$`r_\sigma (t)={\displaystyle \frac{N_\sigma (t)}{N_\sigma (t_d)}},`$
$`r_\pi (t)={\displaystyle \frac{N_\pi (t)}{N_\pi (t\mathrm{})}}=1r_\sigma (t),`$ (23)
and the thermal and nonthermal enhancements $`r_{th}(t)`$ and $`r_{nt}(t)`$ as functions of time scaled by the beginning time $`t_d`$ in the case of $`\mu =0.2GeV`$. The dashed line located at the freeze-out time $`t_f`$ separates the left thermal region and the right nonthermal region. $`48\%`$ of the $`\sigma `$’s decay before the freeze-out, the rest decay into nonthermal pions after the freeze-out. The measurable thermal and nonthermal enhancements compared with the number of direct pions are $`32\%`$ and $`30\%`$, respectively.

The measurable thermal, nonthermal and total pion enhancements

$`r_{th}(t_f)=2{\displaystyle \frac{N_\sigma (t_d)}{N_\pi ^{dir}(t_f)}}(1r_\sigma (t_f)),`$
$`r_{nt}(t\mathrm{})=2{\displaystyle \frac{N_\sigma (t_d)}{N_\pi ^{dir}(t_f)}}r_\sigma (t_f),`$
$`r_{th}(t_f)+r_{nt}(t\mathrm{})=2{\displaystyle \frac{N_\sigma (t_d)}{N_\pi ^{dir}(t_f)}}`$ (24)
are shown in Fig. 5 as functions of chemical potential $`\mu `$. At a first-order chiral transition point any particle undergoes a mass jump. At the low end of this jump, $`m_\sigma <2m_\pi `$, $`\sigma `$ decay is forbidden, while at the high end, $`m_\sigma >2m_\pi `$, the decay begins with a nonzero rate. This sudden decay was shown in Fig.3. On the other hand, the large $`\sigma `$ mass at the beginning of decay means a strong suppression of the sigma number $`N_\sigma (t_d)`$ which dominates the amplitudes of the enhancements, it drops down exponentially with increasing $`\sigma `$-mass. Compared with the decay happen in the crossover region where the decay begins with $`m_\sigma =2m_\pi `$, the thermal and nonthermal pion enhancements are both washed away by the first-order chiral transition. This is clearly displayed in Fig.5. The enhancements start at the endpoint $`E`$ of the first-order transition line and develop in the crossover region. In a wide chemical potential interval, the total pion enhancement is very strong, $`r_{th}(t_f)+r_{nt}(t\mathrm{})0.5`$. The shapes of the enhancements are mainly controlled by the decay rate shown in Fig.3. The nonthermal and total enhancements have their maximum values in the vicinity of the endpoint $`E`$, this is in agreement with our qualitative analysis in the last section.

The momentum distribution of the nonthermal pions is shown in Fig.6 for four different chemical potentials. The nonthermal pions are always centralized in a very narrow momentum region around $`p0.26GeV`$. This narrow distribution leads to a sharp peak in the momentum spectra for the total pions, shown in Fig.7. For low chemical potentials the nonthermal momenta and the most appropriate thermal momenta almost coincide, while for high chemical potentials the nonthermal peak is separated from the thermal peak obviously. The thermal peak is always smooth, but the nonthermal peak is always sharp.

The mean momentum of the nonthermal pions is shown in Fig.8 as a function of chemical potential. Unlike the strong chemical potential dependence of the enhancements, especially around the endpoint $`E`$, the mean momentum is almost a constant, $`\overline{p}_{nt}0.26GeV`$. This value is not so small as the people expected before.

Our study on $`\sigma `$ decay depends on both the chiral dynamics which determines the decay rate, and the expansion mechanism which is reflected in the time evolution of the temperature. To see which one is dominant in the above calculations, we have chosen different expansion parameters. We changed $`\mu _{max}`$ to control the chemical potential dependence of the expansion rate, and the beginning time $`t_d`$ of the decay. While the amplitudes of the measurable thermal and nonthermal enhancements change remarkably, their shapes, the total enhancement and the mean momentum are insensitive to the expansion parameters. The chiral symmetry restoration plays a crucial rule in these observable quantities. We even let the freeze-out temperature $`T_f`$ be the decay temperature $`T_d`$, this means no thermal pions produced by $`\sigma `$ decay, the mean momentum in this case decreases by less than $`10\%`$ only.

## 4 conclusions

Since the $`\sigma `$-mass goes up and the $`\pi `$-mass drops down as the hadronic system cools down, the more and more large mass difference between $`\sigma `$ and $`\pi `$ leads to $`\sigma `$ decay into two pions when the $`\sigma `$-mass exceeds the threshold $`2m_\pi `$. These pions result in thermal enhancement when they are produced before the freeze-out and nonthermal enhancement when they are produced after the freeze-out. Based on the two-flavor NJL model we have presented a detailed analysis of the pion production at finite temperatures and densities, and discussed the relation between the decay and the chiral symmetry restoration, especially the crucial effect of the order of the chiral transition on the thermal and nonthermal enhancements.

We have found that the chiral symmetry breaking and restoration control the pion enhancements of the system. At high densities the number of $`\sigma `$’s at the beginning of their decay is strongly suppressed by the large $`\sigma `$-mass, originating from the first-order chiral transition. Compared with the large number of direct pions produced by the freeze-out of the thermal pion gas, there are almost neither thermal nor nonthermal pion enhancements. In practice, the enhancements begin at the endpoint $`E`$ of the first-order chiral transition line and present only in the crossover region with high temperatures and low densities. The nonthermal enhancement in this region develops very fast in the early stage of the decay and reaches its maximum value in the vicinity of the endpoint $`E`$. For a wide chemical potential interval, $`0\mu 0.25GeV`$, both the thermal and nonthermal enhancements are remarkably large due to the almost constant decay rate during the whole decay process. The maximum enhancement is $`63\%`$ which agrees well with the value $`2/3`$ estimated by simple statistics.

While the total enhancement and the structure of the thermal and nonthermal enhancements seem insensitive to the expansion parameters we have chosen in the calculations, the amplitudes of the thermal and nonthermal enhancements depend strongly on the expansion mechanism. On the other hand, our results are derived in a particular chiral model. Although we expect the characteristics of the enhancements to be of more general validity, the numerical results are of course model dependent. Therefore, more reliable conclusions need more realistic expansion mechanism and more serious QCD model.

acknowledgments: This work was supported in part by the NSFC under grant numbers 19845001 and 19925519, and the Major State Basic Research Development Program under contract number G2000077407.

Figure Captions

Fig.1: The Feynman diagram for $`\sigma `$ decay into $`2\pi `$ to the lowest order in $`1/N_c`$ expansion. The solid lines denote quarks and antiquarks, and the dashed lines denote mesons.

Fig.2: The decay rate $`\mathrm{\Gamma }_{\sigma 2\pi }`$ and the coupling strength $`g_{\sigma \pi \pi }`$ as functions of $`T`$ at $`\mu =0`$. The thick and thin solid lines indicate the decay rates with and without consideration of Bose-Einstein statistics for the final state pions.

Fig.3: The decay rate $`\mathrm{\Gamma }_{\sigma 2\pi }`$ as a function of $`T`$ for different chemical potentials. The lines marked with $`a,b,c,d,e`$ correspond to $`\mu =0,0.1,0.2,0.25,0.3GeV<\mu _E=0.325GeV`$, respectively, and the line $`f`$ is with $`\mu =0.33GeV>\mu _E`$.

Fig.4: The sigma number and pion number scaled by the corresponding maximum value, and the thermal and nonthermal pion enhancements as functions of time scaled by the beginning time of the decay, at $`\mu =0.2GeV`$. The dashed line separates the left thermal and right nonthermal regions.

Fig.5: The measurable thermal, nonthermal and total pion enhancements as functions of $`\mu `$.

Fig.6: The normalized momentum distribution of the nonthermal pions for four different chemical potentials.

Fig.7: The normalized momentum spectra of the total pions for two different chemical potentials.

Fig.8: The mean momentum of the nonthermal pions as a function of $`\mu `$.


format: text

text:
# Pulse interaction in nonlinear vacuum electrodynamics

## 1 Introduction

Although classical electromagnetic theory deals with linear Maxwell equations, there have been numerous attempts to bring the nonlinear phenomena into the stage. All relativistic and gauge invariant versions of electromagnetism are based on the Lagrangian density, $`L`$, which depends on the invariants of the field tensor. Generally, in terms of the electric ($`𝐄`$) and magnetic ($`𝐁`$) fields the Maxwell equations in absence of external charges may be written in a standard form:

$`𝐃_𝐭\times 𝐇=\mathrm{𝟎},𝐃=\mathrm{𝟎},`$
(1)
$`𝐁_𝐭+\times 𝐄=\mathrm{𝟎},𝐁=\mathrm{𝟎},`$
where we put $`c=1`$ and $`𝐃=\frac{𝐋}{𝐄},𝐇=\frac{𝐋}{𝐁}`$. The Lagrangian $`L(I,J^2)`$ depends on Poincare invariants $`I=𝐄^\mathrm{𝟐}𝐁^\mathrm{𝟐}`$ and $`J=\mathrm{𝐄𝐁}`$ only. The distinctive feature of Eqs. (1) is that since the Poincare invariants are identically zero for the plane electromagnetic wave, the latter is insensitive to vacuum nonlinearity and propagates without distortion.

Of particular interest are the nonlinear corrections to the linear electrodynamics arising due to vacuum polarization in the strong electromagnetic field. In the ultimate case of slowly varying fields this results in Heisenberg- Euler electrodynamics, which is discussed in many textbooks (e.g. ).

The main point of this paper is to describe the simplest, in a sense, nonlinear vacuum process: the interaction of two electromagnetic waveforms propagating in opposite directions.

## 2 Maxwell equations

We consider a linearly polarized wave propagating in the $`z`$ direction of the form $`E_x=E(z,t),B_y=B(z,t)`$ with all other components being zero. In this situation, the second Poincare invariant vanishes, $`J0`$, so the Maxwell equations are written as

$`\left(EL(I)_I\right)_t+\left(BL(I)_I\right)_z`$ $`=`$ $`0,`$
$`\left(B\right)_t+\left(E\right)_z`$ $`=`$ $`0,`$
where the subscript denotes the derivative with respect to the corresponding variable and $`I=E^2B^2`$. The Lagrangian in Eq. (LABEL:meq) is expanded in powers of $`I`$. Keeping the lowest-order nonlinear corrections we have

$$L(I)=I+\frac{1}{2}\sigma I^2+\mathrm{}$$
(3)
With the help of the appropriate scale transform, the coefficient $`\sigma `$ may be reduced to $`\pm 1`$. For the particular case of the Heisenberg-Euler electrodynamics, $`\sigma =1`$ . Of interest also is to keep in mind the Born-Infeld electrodynamics (e.g.  ) with the Lagrangian

$$L_{BI}(I)=1\sqrt{1I}.$$
(4)
## 3 Energy-momentum tensor

The conservation laws for Eqs. (LABEL:meq) are given by

$$W_t+N_z=0,N_t+P_z=0,$$
(5)
where the components of the energy-momentum tensor, namely, the energy density, $`W`$, the momentum density, $`N`$, and the stress, $`P`$, may be obtained using standard variation procedure (e.g. ). Explicitely,

$`W`$ $`=`$ $`2E^2L_IL`$
$`N`$ $`=`$ $`2EBL_I`$ (6)
$`P`$ $`=`$ $`2B^2L_I+L.`$
Usually Eqs. (5,6) are thought of as a consequence of the Maxwell equations (LABEL:meq). However, we may consider the relations (6) as a constraint implied upon the components of the momentum-energy tensor, so there are two independent variables in Eqs. (5), for example, $`W`$ and $`N`$. One can easily check that for the nontrivial solutions of Eqs. (LABEL:meq), i.e. for $`I0`$, the Jacobian of the transform $`E,BW,N`$ is non-zero. Thus, instead of looking for the solutions of Eqs. (LABEL:meq) we can solve Eqs. (5,6) excluding the Poincare invariant $`I`$ from Eqs. (6).

## 4 Solution

To exclude $`I`$ it is convenient to introduce the invariants of the energy- momentum tensor, that is, its trace, $`S=PW`$, and the determinant $`T=WPN^2`$. As it follows from Eqs. (6)

$`S`$ $`=`$ $`2(LIL_I)`$
$`T`$ $`=`$ $`I\left(L^2\right)_IL^2.`$
The latter relations implicitly define the dependence $`T=T(S)`$. Substituting the Lagrangian (3) into Eqs. (LABEL:inv) we find that the first nonvanishing term of the expansion of $`T`$ in powers of $`S`$ is linear and it is provided by the quadratic term of the expansion (3): $`T(S)=\sigma S+\mathrm{}`$. It is noteworthy that the Born-Infeld Lagrangian (4) yields exactly the linear dependence $`T(S)=S`$.

The relations (5) are resolved introducing the potential $`\psi `$: $`W=\psi _{zz},N=\psi _{zt},P=\psi _{tt}`$. Restricting ourself with the linear relation between $`T`$ and $`S`$, we obtain the Ampere-Monge type equation for $`\psi `$:

$$\psi _{zz}\psi _{tt}\psi _{zt}^2=\sigma (\psi _{zz}\psi _{tt}).$$
(8)
There are trivial solutions to this equation $`\psi (z,t)=F(z\pm t)`$ with an arbitrary function $`F`$, which correspond to the plane electromagnetic waveforms described by Eqs. (LABEL:meq) with $`I=0`$. Besides these, implementing the Legendre transform  one can easily obtain the general integral of Eq. (8) valid for $`T0`$ and, consequently, for $`I0`$. As a result, we get the components of the energy-momentum tensor in a parametric form:

$`W`$ $`=`$ $`\sigma (F_1^{}(\xi )+F_2^{}(\eta )+2F_1^{}(\xi )F_2^{}(\eta ))/\mathrm{\Delta }(\xi ,\eta ),`$
$`P`$ $`=`$ $`\sigma (F_1^{}(\xi )+F_2^{}(\eta )2F_1^{}(\xi )F_2^{}(\eta ))/\mathrm{\Delta }(\xi ,\eta ),`$
$`N`$ $`=`$ $`\sigma (F_2^{}(\eta )F_1^{}(\xi ))/\mathrm{\Delta }(\xi ,\eta ),`$ (9)
$`z`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\xi +\eta F_1(\xi )F_2(\eta )),`$
$`t`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\xi \eta +F_1(\xi )F_2(\eta )),`$
where $`F_{1,2}`$ are arbitrary functions and $`\mathrm{\Delta }(\xi ,\eta )=1F_1^{}(\xi )F_2^{}(\eta )`$.

Consider, for example, two localized pulses of the arbitrary shape propagating in opposite directions. This corresponds to the following initial conditions:

$$W(z,t)|_t\mathrm{}=W_1(z+t)+W_2(zt),$$
(10)
where $`W_{1,2}(\xi )|_{\xi \pm \mathrm{}}0`$. This initial condition is provided by the following choice of $`F_{1,2}`$ in Eqs. (9):$`F_{1,2}^{}(\xi )=\sigma W_{1,2}(\xi )`$ and $`F_1(\xi )|_\xi \mathrm{}0`$, $`F_2(\xi )|_\xi \mathrm{}0`$. The asymptotic of the solution (9) at $`t\mathrm{}`$ is then given by

$$W(z,t)=W_1(z+t\sigma K_1)+W_2(zt+\sigma K_1),$$
(11)
where

$$K_{1,2}=\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\xi W_{1,2}(\xi )$$

is the net energy carried by the corresponding pulse.

Typical plots $`W(z,t)`$ are depicted at Fig. 1,2. The first picture shows the interaction of two identical pulses. The interaction of one pulse with a sequence of two pulses is shown in Fig. 2.

## 5 Discussion

Of interest is the geometrical sense of the obtained solution (9). The parameters $`\xi `$ and $`\eta `$ are, in fact, the light-cone coordinates disturbed by the electromagnetic field. One may say that the electromagnetic field alters the space-time metric due to the dependence of the speed of light on the field strength. In contrast with general relativity, the space-time remains flat.

Another interesting point is that for $`\sigma =1`$ the increase in the pulse amplitude results in delay in energy (and information) exchange between distant points, that is, the solution described by (11) is subluminal. This takes place for both the Heisenberg-Euler electrodynamics, which is currently the only one of physical sence, and for the elegant Born-Infeld theory, for which our results are exact. However, for $`\sigma =1`$ the pulse propagation would be superluminal.

From the viewpoint of nonlinear physics, the electromagnetic pulses in vacuum exhibit the soliton-like behavior: the collision results in a phase shift but the form of a pulse remains unchanged. The main interesting point with this respect is that unlike usual nonlinear equations, the shape of the soliton is arbitrary.


text:
# Electrically Tunable Collective Modes in a MEMS Resonator Array

## Abstract

Using optical diffraction, we study the mechanical vibrations of an array of micromechanical resonators. Implementing tunable electrostatic coupling between the suspended, doubly-clamped Au beams leads to the formation of a band of collective vibrational modes within these devices. The evolution of these modes with coupling strength is clearly manifested in the optical diffraction pattern of light transmitted through the array. The experimental results are analyzed using a simple model for one-dimensional phonons. These structures offer unique prospects for spectral analysis of complex mechanical stimuli.

Interaction between light and mechanically vibrating systems was first discussed by Einstein and Bohr. A gedanken experiment of optical diffraction by a vibrating two slit structure was employed to demonstrate the complementarity principle of quantum mechanics . These ideas were later elaborated in order to formulate the theory of x-rays and neutrons diffraction by crystals  . In these systems mechanical vibrations due to both thermal and quantum fluctuations strongly modify the diffraction pattern. This allows experimental study of the crystal’s normal modes of vibration (phonons) and their dispersion relation. Micro electro-mechanical systems (MEMS) technology allows studying similar phenomena using artificial mechanical systems at mesoscopic length scales. Such studies are motivated not only by scientific interest but rather also by the prospects of developing new micro opto-mechanical devices. This rapidly growing field of research employs MEMS technology to realize a variety of on-chip fully integrated optical devices (for a review see ). A particular example related to the present work is the movable diffraction grating. In this device the efficiency of diffraction is controlled by moving the grating with respect to the substrate beneath it, allowing thus optical modulation   . Operation of such devices at relatively high frequencies however requires a good understanding of the dynamics and mechanical properties of the system.

In the present paper we use optical diffraction to study the mechanical properties of a periodic array (grating) of suspended doubly-clamped beams made of Au. What is novel and especially interesting about the present work is that our devices allow application of mutual electrostatic forces between the beams. This coupling  gives rise to the formation of a band of collective modes of vibration (phonons). We excite these collective modes parametrically and employ optical diffraction to study the response. A simple model describing our system is presented and compared with experiment. We conclude by briefly describing some unique device applications of this system.

The bulk micromachining process employed for sample fabrication is described in Fig. 1. In this process the substrate beneath the grating is completely etched away, thus allowing optical access to the grating from both sides. In the first step chemical vapor deposition is employed to deposit a 70 $`nm`$ thick layer of low-stress silicon nitride on the front and back sides of a Si wafer. A square window is opened in the silicon nitride on the back side using photolithography and wet etching (Fig. 1(a)). The structure shown in Fig. 1(b) is realized through a highly selective, anisotropic KOH etch for the backside of the Si wafer. This occurs within the patterned region and yields a 270$`\mu `$m square silicon nitride membrane on the front side of the wafer. The grating beams and adjacent electrodes are fabricated on top of this membrane using electron beam lithography, followed by thermal evaporation of Au and liftoff (Fig. 1(c)). Each resulting beam has length $`b=270`$ $`\mu `$m, width $`1\mu `$m and thickness $`0.25\mu `$m (measured using an atomic force microscope) and the grating period is $`a=4\mu `$m. In the final step the membrane is removed using electron cyclotron resonance (ECR) plasma etching from the back side of the sample. This process step employs an Ar/NF3 gas mixture, and results in suspension of the Au beam array (Fig. 1(d)). Figure 1(e) shows a side view micrograph of the device. The electrodes form two interdigitated combs; with fingers alternately connected to the two base electrodes. This design allows application of electrostatic forces between the beams.

To characterize uniformity within the device we measure the fundamental resonance frequency of each suspended beam in the array. This is done in-situ, using the output from a commercial scanning electron microscope’s imaging system to detect mechanical displacement. We have employed this technique previously to study the mechanical properties of individual, similarly-fabricated Au beams  . Fig. 2(a) shows a typical response peak from an individual beam. By measuring all 67 suspended beams in the array we find that the distribution of resonance frequencies has an average of 179.3 kHz and a standard deviation of 0.53 kHz (see Fig. 2(b)). Individual mechanical quality factors $`Q`$ range from 2,000 to 10,000. Note that no correlation is found between the location of the beam within the device and its specific resonance frequency or $`Q`$; the small beam-to-beam variations appear to be random.

What is expected when a voltage $`V`$ is applied between the two combs? We employ a simple one-dimensional model for an $`N`$-element array of coupled pendulums  to describe our system (see Fig. 3(a)). While the first and last pendulums in the array are clamped and stationary, all others ($`n=2,3,\mathrm{},N1`$) are free to oscillate about their equilibrium positions $`na`$. Here $`a`$ represents the equilibrium spacing between neighboring pendulums. In the absence of any coupling, the angular frequency for small oscillations of each (identical) pendulum is $`\omega _0`$. The displacement of the system is described by a set of coordinates $`x_n`$ ($`n=1,2,\mathrm{},N`$) (see Fig. 3(a)). Applying a voltage $`V`$ gives rise to an attractive interaction between each pendulum and its nearest neighbors $`\phi \left(s\right)=C\left(s\right)V^2/2`$, where $`s`$ is the distance between the interacting pendulums and $`C`$ is the capacitance. Neglecting coupling between non-neighboring pendulums, and assuming small oscillations, we find the following set of equations of motion:

$$m\underset{n}{\overset{}{x}}=m\omega _0^2x_n+u(2x_nx_{n1}x_{n+1}),$$
(1)
where $`n=2,3,\mathrm{},N1`$, and $`u=\phi ^{^{\prime \prime }}\left(a\right)`$ (here dots represents time derivatives and primes represent spatial derivatives). Note that $`u>0`$ due to the attractive nature of the interaction between nearest neighbors. These equations can be greatly simplified by employing a transformation to the eigenmodes (phonons) of the system:

$$x_n=\underset{m=2}{\overset{N1}{}}\chi _n^{(m)}v_m,$$
(3)
where $`\chi _n^{(m)}=\sqrt{2/\left(N1\right)}\mathrm{sin}\left(k_m\left(n1\right)a\right)`$ is the spatial shape of mode number $`m`$ ($`m=2,3,\mathrm{},N1`$), $`k_m=\left(m1\right)\pi /L`$ is the wavevector, and $`L=\left(N1\right)a`$ is the length of the system. Figure 3(b) shows the shape of the three lowest modes $`m=2,3,4`$. Substituting Eq. (3) in Eq. (1) leads to a set of decoupled equations of motion:

$$\underset{m}{\overset{}{v}}=\omega _m^2v_m,$$
(4)
where $`\omega _m^2\left(V\right)=\omega _0^2\left(2C^{^{\prime \prime }}\left(a\right)/m\right)V^2\mathrm{sin}^2\left(k_ma/2\right)`$. A stationary voltage $`V=V_{dc}`$ thus gives rise to the formation of a band of collective modes between frequencies $`\omega _0`$ and $`\omega _b=\sqrt{\omega _0^2\left(2C^{^{\prime \prime }}\left(a\right)/m\right)V_{dc}^2}`$. The associated wavevectors, $`k_m`$, vary from zero to $`\pi /a`$ (see Figure 3(c)).

Each mode can be selectively excited by adding an AC voltage to the DC bias, namely $`V=V_{dc}+V_{ac}\mathrm{cos}\left(\gamma t\right)`$. Assuming $`V_{ac}<<V_{dc}`$ we find from Eq. (4):

$$\underset{m}{\overset{}{v}}=\omega _m^2\left(V_{dc}\right)(1h_m\mathrm{cos}\left(\gamma t\right))v_m,$$
(5)
where $`h_m=2\left(V_{ac}/V_{dc}\right)\left[\left(\omega _0/\omega _m\right)^21\right]`$. Thus an AC voltage component gives rise to parametric excitation of each mode with amplitude $`h_m`$   . Parametric resonance occurs when the frequency of the AC voltage $`\gamma `$ is close to $`2\omega _m/l`$, where $`l`$ is an integer. Near these values the system may exhibit unstable behavior in which the amplitude of oscillations grows as a function of time. In the linear theory of parametric resonance this growth is exponential, and occurs when the amplitude of parametric excitation $`h`$ exceeds a critical value dependent upon the damping in the system. Most systems, however, possess some degree of non-linearity which comes into play as soon as the amplitude of the motion becomes appreciable. Thus, while the linear theory is useful in determining the conditions for the occurrence of parametric resonance, it is inadequate for determining the steady-state response of the system. Unfortunately, the nonlinear coefficients of our system are not known and therefore its steady-state response cannot be predicted. However, we expect that the response of modes with high index $`m`$ will be relatively large because $`h_m`$ increases in magnitude with $`m`$.

We detect the collective mechanical vibrations of the array by diffraction measurements. With spatially-uniform light incident upon the array, the intensity of diffraction is proportional to the form factor $`\left|\varphi \left(q_x\right)\right|^2`$ where:

$$\varphi \left(q_x\right)=\underset{n=1}{\overset{N}{}}\mathrm{exp}\left[iq_x\left(na+x_n\left(t\right)\right)\right].$$
(6)
Here $`q_x`$ is the $`x`$ component of the change in wavevector between incoming and outgoing waves (the $`x`$ direction lies in the plane of the sample and is perpendicular to the long axis of the beams) . Consider the case where the system is tuned to a diffraction peak, namely $`q_xa=2\pi l`$ with $`l`$ an integer. As shown below, for this particular case the interpretation of the experimental results becomes greatly simplified. We calculate the form factor for the case where mode $`m`$ oscillates with amplitude $`A_m`$ and all other modes are stationary. Assuming small oscillations, namely $`q_xx<<1`$, we find:

$$\left|\varphi \left(q_x\right)\right|^2=N^2+\left[q_xA_m\mathrm{cot}\left(k_ma/2\right)\mathrm{cos}\left(\omega _mt\right)\right]^2\frac{\left(1\right)^m+1}{N1}.$$
(7)
Thus the collective mechanical oscillations give rise to an oscillatory component in the diffraction signal at angular frequency $`2\omega _m`$. Note that with uniform illumination any mode with $`m`$ odd will not contribute to diffraction. This occurs due to the exact cancelation of the response from different portions of the entire device. Similar cancelation gives rise to the $`\mathrm{cot}^2\left(k_ma/2\right)`$ factor in Eq. (7).

The optical setup utilized for the diffraction measurements is schematically depicted in Fig. 4. A polarization maintaining single mode fiber (numerical aperture NA$`=0.15`$ and core diameter $`d=9`$ $`\mu `$m) delivers infrared light from a laser to the sample. We use a tunable wavelength diode laser operating in the range 1535-1635 nm. A spherical lens (focal length $`f=1`$ mm) collimates the beam, which illuminates the back side of the grating at normal incidence with respect to the array plane. At the array plane, located a distance $`s_1=3`$ cm from the collimator, we estimate the diameter of the beam to be $`2f`$NA$`+s_1d/f=570`$ $`\mu `$m. For the results presented, the polarization of the incident light is TE, i.e. the electric field is parallel to the grating lines. We note, however, that similar results (not presented here) were obtained with perpendicular polarization. For the present case there are four diffraction peaks with angles $`\theta _n=\mathrm{sin}^1\left(n\lambda /a\right)`$ ($`n=\pm 1,\pm 2`$) with respect to the normal incidence. Here $`\lambda `$ is the wavelength of incident light.. Transmitted light is collected by a cylindrical lens and focused into a second, single-mode fiber. The lens enables the intensity of collected light to be maximized without degrading the spectral resolution. The fiber delivers light collected from the first order ($`n=1)`$ diffraction peak to a photodiode detector. The distance between the sample and the lens, $`s_2=1`$ cm, was chosen as a compromise between two conflicting considerations, namely, simultaneously maximizing spectral resolution and light intensity. While the former consideration favors a large distance, the later favors a small one. For the distance chosen, the spectral resolution is $`\delta \lambda =a\delta \theta \mathrm{cos}\theta _1=a\left(d/s_2\right)\mathrm{cos}\theta _13`$ nm. All measurements are done at room temperature in vacuum of $`10^3`$ torr.

We first study diffraction from the array in the absence of any interelectrode bias voltage. The inset of Fig. 4 shows the intensity detected by the photo diode as a function of $`\lambda `$. The full width at half maximum (FWHM) of the diffraction peak is estimated using the Fraunhofer diffraction formula  to be $`\delta \lambda \lambda /N=23`$ nm. We find good agreement between this estimate and the measured value. No effect is resolved from the thermally-driven mechanical vibrations of the beams. Their effect upon diffraction can be characterized by multiplying the diffracted intensity by a Debye-Waller factor, $`\mathrm{exp}\left(2W\right)`$  . In the present case we estimate that $`2W10^8`$, hence thermal fluctuations are not expected to affect the diffraction from our micromechanical array significantly. This is confirmed by the experimental data.

In order to excite the modes of vibrations of the system externally, we apply a voltage $`V=V_{dc}+V_{ac}\mathrm{cos}\left(\gamma t\right)`$ between either side of the interdigitated electrode arrays. We tune the laser wavelength to the diffraction peak ($`\lambda =1582`$ nm) and measure the photodetector response using a lock-in amplifier operating in $`2f`$ mode. This allows us to detect the Fourier component of the array response at angular frequency $`2\gamma `$. Figure 5 is a color map plot showing this second harmonic response, $`R,`$ as a function of both $`V_{dc}`$ and $`f=\gamma /2\pi `$. The amplitude of the AC voltage is $`V_{ac}=50`$ mV for these measurements.

With $`V_{dc}=0`$ we find a peak in $`R`$ at $`f=179.3`$ kHz, associated with the fundamental frequency of the decoupled beams. The FWHM of the peak is 0.6 kHz, close to the standard deviation found in the distribution of the measured fundamental frequencies. This leads us to the conclusion that the width of the response peak at $`V_{dc}=0`$ is dominated by inhomogeneous broadening caused by the non-uniformity of the array.

As we increase $`V_{dc}`$ we observe a gradual increase in the frequency range where relatively large response is observed; we associate this with the formation of a band of collective modes. The lower frequency bound of this range $`f_b=\omega _b/2\pi `$ (the bottom of the band) for relatively small $`V_{dc}`$ is given theoretically by $`f_b=\left(\omega _0/2\pi \right)\left(1C^{^{\prime \prime }}\left(a\right)V_{dc}^2/m\omega _0^2\right)`$. A least squares fit to the measured data (see dashed line in Fig. 5) yields $`C^{^{\prime \prime }}\left(a\right)/m\omega _0^2=2.7\times 10^4`$ V-2. For comparison we derive below a rough order of magnitude estimate of this factor. We substitute $`m`$, which should represent an effective mass of each beam, by the actual mass and we use the approximation $`C^{^{\prime \prime }}\left(a\right)\epsilon _0b/g^2`$, where $`\epsilon _0`$ is permittivity of free space and $`g`$ is the gap between neighboring beams. These crude approximations yield a value of $`1.6\times 10^4`$ V-2 which is quite close to the value deduced from the experimental data. The upper frequency bound of the band, on the other hand, shows some discrepancy with theory. While the measured value depends on $`V_{dc}`$, our simple model predicts a upper value that is fixed. We obtained similar behavior of this upper frequency with orthogonal polarization, and with another sample of similar design.

The rich and detailed structure of the frequency-dependent response $`R`$ observed in our experiments is not fully understood. Experimentally, Figure 5 shows that $`R`$ is oscillatory as a function of frequency, and has a relatively large magnitude close to the lower limit of the band. Theoretically this frequency dependence can be found from Eq. (7) and the density of states (which is determined by the dispersion relation). This leads to an approximated proportionality relation $`RA_m^2\mathrm{cos}\left(k_ma/2\right)/\mathrm{sin}^3\left(k_ma/2\right)`$ for even $`m`$. As discussed above, however, our simple linear theory is not adequate for predicting the steady-state response $`A_m`$. Comparison with experiment leads to the following conclusions: (a) For small $`m`$ $`A_m`$ goes to zero faster than $`m^{3/2}`$. (b) Near the lower bound of the band ($`k_ma\pi `$) $`A_m`$ grows faster than $`\sqrt{1/\mathrm{cos}\left(k_ma/2\right)}`$. (c) The spiky behavior indicates that $`A_m`$ is not a monotonic function of $`m`$. These intriguing details of the device response warrant further investigation.

Electrically tunable arrays offer unique prospects for opto-mechanical signal processing devices such as tunable filters and optical modulators. A particularly intriguing example is opto-mechanical spectral analysis of electrical waveforms. Consider an arbitrary electrical signal applied between the two comb electrodes of our device. Its Fourier components falling within the vibrational band (formed by the DC electrostatic coupling as described above) will parametrically drive the collective modes of the array. Each of these excited modes will result in a diffracted order with strength directly proportional to the respective Fourier component. Since each order is diffracted at a characteristic angle, continuous, real-time spectral analysis of the applied waveform can be realized using a photodetector array. In principle, by scaling the size of the resonant beams downward into the realm of NEMS (nano electromechanical systems), the operating frequency of such a device can be extended to very high frequencies .

In conclusion, we have demonstrated the ability to induce and control collective modes of mechanical vibration within an artificial mesoscopic lattice. Further experimental and theoretical work will elucidate what appear to be rich dynamics in such systems, and will allow their optimization for novel micro opto-mechanical device applications.

The authors are grateful to K. Schwab for his assistance in sample fabrication. This research was supported by DARPA MTO/MEMS under grant DABT63-98-1-0012. E. B. gratefully acknowledges support from the Rothschild fellowship and the R. A. Millikan postdoctoral fellowship at Caltech.