text:
# Transverse Momentum Spectra of 𝐽/𝜓 and 𝜓' Mesons from Quark Gluon Plasma Hadronization in Nuclear Collisions

## Abstract

Recent results on transverse mass spectra of $`J/\psi `$ and $`\psi ^{}`$ mesons in central Pb+Pb collisions at 158 A$``$GeV are considered. It is shown that those results support a hypothesis of statistical production of charmonia at hadronization and suggest the early thermal freeze–out of $`J/\psi `$ and $`\psi ^{}`$ mesons. Based on this approach the collective transverse velocity of hadronizing quark gluon plasma is estimated to be $`v_T^H0.2`$. Predictions for transverse mass spectra of hidden and open charm mesons at SPS and RHIC are discussed.

Key words: statistical production of charmonia, early freeze–out, transverse flow, inverse slope parameter

The idealized concepts of chemical (hadron multiplicities) and thermal (hadron momentum distributions) freeze-outs were introduced to interpret data on hadron production in relativistic nucleus–nucleus (A+A) collisions . The first experimental results on yields and transverse mass ($`m_T=\sqrt{p_T^2+m^2}`$) spectra suggested the following scenario: for the most abundant hadron species ($`\pi ,N,K,\mathrm{\Lambda }`$) the chemical freeze-out, which seems to coincide with the hadronization of the quark gluon plasma (QGP), is followed by the thermal freeze–out occurring at a rather late stage of the A+A reaction. In this letter we discuss whether the new data of NA50  on transverse mass spectra of $`J/\psi `$ and $`\psi ^{}`$ mesons produced in central Pb+Pb collisions at 158 A$``$GeV are consistent with the above picture. Our consideration is based on a recent hypothesis  of statistical production of charmonia at hadronization. We further suggest that thermal freeze–out of charmonia coincides with the hadronization transition. The consequences of this assumptions are in agreement with existing data on $`m_T`$ spectra. They allow to extract the transverse flow velocity of the hadronizing QGP and lead to new predictions which can be tested by future measurements.

Rescattering among partons and, in the late stage of the reaction process, hadrons created in relativistic A+A collisions should cause local thermalization and the development of transverse collective flow of matter. Thus the final transverse motion of hadrons can be considered as a convolution of the transverse flow velocity of the freezing–out matter element with the thermal motion of the hadrons in the rest frame of this element. The resulting $`m_T`$ spectrum has approximately exponential shape:

$$\frac{1}{m_T}\frac{dn}{dm_T}Ce^{m_T/T^{}},$$
(1)
where $`C`$ and $`T^{}`$ are a normalization and an inverse slope parameters, respectively. The $`T^{}`$ parameter is related to the thermal freeze–out temperature $`T_f`$ and the mean transverse flow velocity $`v_T^f`$. In the nonrelativistic approximation ($`m_T<2m`$) contributions from thermal and collective particle motions can be separated and the inverse slope parameter of the observed hadron spectrum can be expressed as :

$$T^{}=T_f+\frac{2}{\pi }mv_T^f^2.$$
(2)
Note that factor $`2/\pi `$ appears in Eq.(2) due to the cylindrical geometry expected in high energy A+A collisions (see  for details). In the case of the simultaneous thermal freeze–out of all hadrons the inverse slope parameter $`T^{}`$ should follow the linear dependence on $`m`$ given by Eq. (2). In Fig. 1 we present a compilation of inverse slope parameters measured for various hadron species in central Pb+Pb collisions at 158 A$``$GeV. It is observed that the simplistic hypothesis of common thermal freeze–out of all hadrons is not supported by the data.

The freeze–out condition for pions can be determined independently of the other hadron species from the results on two pion correlations. This procedure leads to the following values of the pion thermal freeze–out parameters :

$$T_f=120\pm 12\mathrm{MeV},v_T^f=0.5\pm 0.12.$$
(3)
In Fig. 1 the upper solid line indicates the dependence given by Eq. (2) for the pion freeze–out parameters. The values of $`T^{}`$ for $`K^+`$ and $`p`$ are close to those given by this line (2,3). The results on $`T^{}`$ for other hadrons are in general below the ’pion freeze-out line’. This may be

interpreted in the following way. For hadrons with “low” masses and large interaction cross sections the thermal freeze–out takes place at the late stage of the expansion, i.e., at large $`v_T^f`$ and small $`T_f`$ given by Eq. (3). Heavy hadrons (due to large mass and smaller cross section) are expected to decouple from the system early thus leading to the smaller values of $`T^{}`$ than those expected from the ’pion freeze–out line’ ***For pions, in contrast to other hadrons, the relativistic effects neglected in Eq. (2) are important. This is a main reason why the measured pion inverse slope parameter is above the thermal pion freeze–out line in Fig. 1.. This is seen from the $`m_T`$-slopes of the hyperons and most clearly can be illustrated by a small $`T^{}`$ value of the $`\mathrm{\Omega }`$ hyperon shown in Fig. 1. This explanation is in the line with the typical hydrodynamical calculations (see, for instance, ) – the inverse slope parameter $`T^{}`$ of the hadron is smaller, if it decouples early, i.e., at higher temperatures $`T_f`$, but smaller velocity $`v_T^f`$. A more careful analysis requires the combination of the hydrodynamic approach at the early stage of the expansion with cascade model calculations at the latest stage .

The yields of produced light hadrons are surprisingly well reproduced within the statistical approach to particle production. The temperature parameter extracted from the fit to the multiplicity data is found to be $`T_H=175\pm `$10 MeV . This chemical freeze–out temperature appears to be quite close to the expected temperature of the transition between the hadron gas and the quark gluon plasma. This fact suggests the possibility to ascribe the observed statistical properties of hadron yield systematics at high energies to the statistical nature of the hadronization process, i.e., the hadrons are ”born” at hadronization temperature $`T_H`$ into the state of chemical equilibrium and their multiplicities are frozen–out afterwards.

Within this approach one can make a rough estimate of a lower limit for the value of the measured inverse slope parameters:

$$T_{MIN}^{}=T_H175\mathrm{MeV}.$$
(4)
This is done under two extreme assumptions: there is no collective transverse flow of hadronizing matter and there is no rescattering between produced hadrons. The measured values of $`T^{}`$ parameter for all hadrons satisfy the condition $`T^{}>T_{MIN}^{}`$ showing that the approach sketched above leads to self-consistent results in the light hadron sector.

A very different approach is traditionally used in modeling the production process of heavy hadrons like quarkonia. It is based on the assumption that quarkonia are created significantly prior the hadronization of QGP. This in general implies that the systematics of quarkonium production should be different than the one established for light hadrons. However recently it was found that multiplicity of the $`J/\psi `$ mesons in nuclear collisions at high energies follows dependences well known from light hadron data  and that it is consistent with the hypothesis of statistical production of $`J/\psi `$ mesons at hadronization . In particular the data on $`J/\psi `$ and $`\psi ^{}`$ yields in central Pb+Pb collisions at 158 A$``$GeV are consistent with the predictions of statistical models  for a typical values of $`T_H175`$ MeV extracted from light hadron systematics.

The new hypothesis of statistical $`J/\psi `$ production at hadronization can be further tested using data on $`m_T`$ spectra. As it follows from the above discussion one expects for $`J/\psi `$ mesons an exponential shape of the $`m_T`$ distribution with the inverse slope parameter $`T^{}(J/\psi )>T_H`$. Recently published experimental results of NA50  on transverse mass spectra of $`J/\psi `$ mesons in central Pb+Pb collisions at 158 A$``$GeV confirm this expectation: the shape of the spectrum is approximately exponential with the fitted inverse slope parameter $`T^{}(J/\psi )=245\pm 5`$ MeV. The measured $`T^{}(J/\psi )`$ value is significantly smaller than that expected on the base of the pion freeze–out line ($`T^{}`$ 610 MeV).

The ‘low’ value of $`T^{}`$ for $`J/\psi `$ suggests its rather early thermal freeze–out. One may argue that this is due to the large mass and low interaction cross section of the $`J/\psi `$ meson. We postulate therefore that the thermal freeze–out of $`J/\psi `$ coincides with the hadronization of QGP, i.e., that the $`J/\psi `$ meson does not participate in the hadronic rescattering after hadronization. It is, however, natural to expect that there is a significant collective transverse flow of hadronizing QGP developed at the early stage of partonic rescattering. Consequently the inverse slope parameter of $`J/\psi `$ meson as well as all other hadrons for which chemical and thermal freeze–outs coincide with hadronization can be expressed as:

$$T_H^{}=T_H+\frac{2}{\pi }mv_T^H^2,$$
(5)
where $`v_T^H`$ is the mean transverse flow velocity of the QGP at the hadronization. Assuming $`T_H=175`$ MeV and using measured value of $`T^{}(J/\psi )=245`$ MeV we find from Eq.(5): $`v_T^H0.19`$. As expected the obtained transverse flow velocity of QGP at hadronization is significantly smaller than the transverse flow velocity of pions ($`0.5`$). The linear $`m`$–dependence of $`T_H^{}`$ (5) is shown in Fig. 1 by the lower solid line. Within the approach discussed here, Eq. (5) can be used to obtain a next estimate of the lower limit of the measured inverse slope parameters for all hadrons. In fact the values of the parameter $`T^{}`$ for all light hadrons are higher than $`T_H^{}(m)`$. The recent results  on the $`m_T`$ spectra of the $`\psi ^{}`$ meson indicate that $`T^{}(\psi ^{})T_H^{}(\psi ^{})=258\pm 5`$ MeV, which suggests that also the $`\psi ^{}`$ meson (like $`J/\psi `$) does not participate in the hadronic rescattering.

One may expect that the thermal freeze–out may coincide with hadronization also for the $`D`$ meson. Under this assumption we calculate the value of the apparent temperature for the $`D`$ meson: $`T^{}(D)T_H^{}(D)217`$ MeV. Note that this result is significantly lower than the predictions of the model of Ref.  although in  the same value of the hadronization temperature was used.

Another important issue is the production of open and hidden charm particles in A+A collisions at RHIC energies. One can expect stronger transverse collective flow effects than at the SPS. This will lead to a linear mass dependence (5) of the apparent inverse slope with approximately the same value of $`T_H175`$ MeV, but with a larger value of $`v_T^H`$. It is necessary to mention that a recent analysis  of the particle number ratios of the RHIC data leads to the above value $`T_H175`$ MeV of the chemical freeze–out temperature. The preliminary RHIC data of hadron $`m_T`$ spectra measured by STAR  and PHENIX  collaborations support the fact of a larger value of the transverse velocity. However, the measurements were made for $`\pi ^\pm ,K^\pm ,\overline{p},`$ and $`p`$ only, and, therefore, this corresponds to the late thermal freeze-out stage.

As long as the RHIC data on the inverse slope parameters for the open and hidden charm mesons are absent it is interesting to compare the different model predictions. Assuming $`T_H175`$ MeV we need the value of $`v_T^H`$ to estimate the inverse slope parameters for charmed hadrons in our approach. The hydrodynamic calculations of Ref.  predict the value of $`v_T^H0.30`$ at the hadronization in Au+Au collisions at RHIC. This leads to an increase of the inverse slopes of charmed hadrons at RHIC in comparison to those values at SPS, e.g., $`T^{}(J/\psi )T_H^{}(J/\psi )350`$ MeV and $`T^{}(D)T_H^{}(D)280`$ MeV.

Note that a significantly larger value of $`T^{}(D)380`$ MeV in Au+Au collisions at RHIC would be obtained, if the thermal freeze-out of D mesons happens at temperature of $`130`$ MeV. On the other hand, the calculations done with the PYTHIA generator for $`p+p`$ collisions show an even higher value of the inverse slope for D mesons $`T^{}(D)440`$ MeV . In fact, the situation is even more uncertain as the dynamical transport calculations of the HSD model (Hadron String Dynamics) made for RHIC energy  predict a rather small inverse slope of about 225 MeV for all particle species. Therefore, the $`T^{}(m)`$ systematics at RHIC for $`m_Tm1`$ GeV will provide an opportunity to find whether the statistical hadronization works, whether a thermal freeze-out of charmed particles happens simultaneously with hadronization and whether the transverse collective flow at hadronization of the QGP is stronger at RHIC than at SPS.

In summary, recent results on transverse mass spectra of $`J/\psi `$ and $`\psi ^{}`$ mesons in central Pb+Pb collisions at 158 A$``$GeV are considered. It is shown that these data support the hypothesis of the statistical production of charmonia at hadronization and suggest a simultaneous hadronization and the thermal freeze–out for $`J/\psi `$ and $`\psi ^{}`$ mesons. Based on this approach the collective transverse velocity of hadronizing quark gluon plasma is estimated to be $`v_T^h0.2`$. Prediction for transverse mass spectra of hidden and open charm mesons at SPS and RHIC are discussed.

Acknowledgments

The authors are thankful to S. V. Akkelin, A. Dumitru, R. Pisarski, J. Schaffner-Bielich and Yu. Sinyukov for stimulating discussions. The very fruitful and stimulating discussions with L.D. McLerran are appreciated. K.A.B. gratefully acknowledges the warm hospitality of the BNL Nuclear Theory Group, where parts of this work were done. The financial support of DAAD, Germany, is acknowledged. The research described in this publication was made possible in part by Award No. UP1-2119 of the U.S. Civilian Research & Development Foundation for the Independent States of the Former Soviet Union (CRDF). This manuscript has been authorized under Contract No. DE-AC02-98CH10886 with the U.S. Department of Energy.


format: text

text:
# Equilibrium and non-equilibrium electron tunneling via discrete quantum states

## Abstract

Tunneling is measured via the quantum levels of a metal nanoparticle. We analyze quantitatively the resonance energies, widths, and amplitudes, both in the regime where only one state is accessible for tunneling and in the non-equilibrium regime when additional states are made accessible one-by-one. For tunneling through one state, our results agree with expectations for sequential tunneling, but in the non-equilibrium regime the resonances are broadened and shifted in ways that require taking into account electron interactions and relaxation.

In nanometer scale devices, electron transport can occur through well-resolved quantum states. Transport in semiconductor quantum dots, metal nanoparticles, and molecules can all be understood within a similar framework, in terms of the energies of the states and the rates for transitions between states. Here we make a quantitative investigation of the fundamental processes at work when electrons tunnel through a discrete-state system, using an Al nanoparticle in a single-electron transistor. By varying gate and source-drain voltages ($`V_g`$,$`V`$), we can manipulate electron flow controllably through one state, or through many, and we can extract tunneling rates for each level. When only one state participates in tunneling, the transport properties are in accord with expectations for simple sequential tunneling. However, for larger voltages, the resonance energies, widths, and currents can all be modified by the population of excited quantum states. For an understanding of the high-$`V`$ regime, non-equilibrium transitions beyond those considered previously must be taken into account.

A cross-sectional device schematic is shown in Fig. 1(a). The use of an aluminum particle with aluminum oxide tunnel barriers provides mechanical and charge stability, and allows $`V`$ and $`V_g`$ to be varied without significantly altering barrier resistances. Fabrication is done using electron-beam lithography and reactive ion etching to create a bowl-shaped hole in a silicon-nitride membrane, with a minimum diameter $`10`$ nm. A gate electrode is formed by depositing 18.5 nm Al, followed by anodization to 3.5 V in an oxygen plasma, and then deposition of 8.5 nm of SiOx . The rest of the device is made by depositing a thick Al electrode onto the bowl-shaped side of the membrane, oxidizing for 3 min in 50 mTorr of O2, depositing 1.5 nm Al onto the other side of the device to make a layer of Al nanoparticles, oxidizing, and then depositing the lower Al electrode. Device parameters are determined from the large-$`V`$ structure of the 4.2-K Coulomb staircase curve . The capacitance of nanoparticle to the top electrode is $`C_L=7.9\mathrm{aF}`$, to the bottom electrode $`C_R=2.7\mathrm{aF}`$, and the gate capacitance is $`C_g=0.06\mathrm{aF}`$. The sum of the resistances of

the two tunnel junctions is $`R_\mathrm{\Sigma }3M\mathrm{\Omega }`$, with individual resistances sufficiently large that the intrinsic widths of the quantum states are smaller than $`k_BT`$. Assuming a roughly hemispherical  particle shape, and a capacitance per unit area of $`0.05\mathrm{aF}/\mathrm{nm}^2`$ , we estimate a nanoparticle diameter $`10`$ nm.

In Fig. 1(b), we plot the differential conductance $`dI/dV`$ as a function of $`V_g`$ and $`V`$, when the sample is cooled in a dilution refrigerator with copper-powder filters on the electrical leads. The lines in the figure are due to tunneling resonances through discrete quantum states in the nanoparticle. Lines having a positive slope correspond to tunneling thresholds across the lower-resistance junction L, and negative slopes are thresholds across junction R . The discontinuity evident in the figure is due to a $`V_g`$-driven change in the charge on another nanoparticle adjacent to the one through which tunneling occurs. This merely shifts the electrostatic potential of the current-carrying particle. The intrinsic energies, current levels, and widths of the resonances are not otherwise altered, so that the full $`dI/dV`$-spectrum can be constructed. From the absence of spin-Zeeman splitting in a magnetic field (not shown) for resonance lines I and II, we can identify these transitions with tunneling from an odd number of electrons $`n_0`$ on the particle to an even number . Since these resonances require increased $`|V|`$ as a function of $`V_g`$, they are $`n_0n_01`$ transitions. (See Fig. 2.) Resonances III and IV correspond to even $`(n_01)`$ odd $`n_0`$ transitions. The large gaps in $`V`$ between each of resonances I and II and the next parallel lines are due to the energy difference $`2\mathrm{\Delta }`$ between a fully paired superconducting state in the Al particle and the next lowest-energy tunneling state with 2 quasiparticles .

Equilibrium Regime: We first consider the region of $`V_g`$ near $`650`$ mV, where $`V_g`$ and $`V`$ can be adjusted so that a single spin-degenerate quantum state (state 0) is accessible for tunneling ($`I`$-$`V`$ curves are shown in Fig. 1(c)). For the case under consideration, in which the quantum level is either empty or singly occupied, the current predicted for sequential tunneling  is

$$I=e\frac{2\gamma _{0L}\gamma _{0R}(f_Lf_R)}{(1+f_L)\gamma _{0L}+(1+f_R)\gamma _{0R}}$$
(1)
where $`\gamma _{0L}`$ ($`\gamma _{0R}`$) is the bare rate for an electron to tunnel from the quantum state 0 to an unoccupied density of states in electrode L (R) and $`f_i`$ is the occupation probability for states in electrode $`i`$ with energy equal to the resonance state (for a thermal distribution $`f_i=[1+exp[(ϵ_0\mu _i)/kT]]^1`$, with $`ϵ_0`$ the energy to occupy the quantum state, and $`\mu _i`$ the chemical potential in electrode $`i`$). Our observations are in excellent accord with this model. For instance, the tunneling current through the quantum state is not the same for both bias directions, being $`I_+=16.6\pm 0.1\mathrm{p}A`$ for $`V>0`$ (i.e. crossing lines II or III) and $`I_{}=8.4\pm 0.1\mathrm{p}A`$ for $`V<0`$ (crossing lines I or IV). This has been observed previously  and is a consequence of spin-degeneracy. For $`V>0`$, electrons tunnel across the high-resistance rate-limiting tunnel barrier into an empty state, so that either spin-up or spin-down electrons can tunnel. For $`V<0`$, the rate-limiting step is for an electron of a given spin on the particle to tunnel through the high-resistance barrier, and the current level is approximately cut in half. By equating the measured currents to Eq. (1), we determine $`\gamma _{0R}=(5.3\pm 0.1)\times 10^7s^1`$ for the high-resistance junction

and $`\gamma _{0L}7\times 10^9s^1`$. In addition, for resonances (e.g., III, IV) in which the spin-degeneracy of the state is split by a magnetic field, the currents through the two Zeeman states for a given bias direction are not equal , in agreement with the simple tunneling theory. The maximum current through the lower-energy Zeeman state is $`e\gamma _{0L}\gamma _{0R}/(\gamma _{0L}+\gamma _{0R})`$ for either bias direction, and the second state then adds current to produce the maximum allowed by Eq. (1).

An interesting feature of Eq. (1) when the tunneling threshold is across the lower resistance barrier (peaks II and IV) is that the maximum of $`dI/dV`$ does not occur exactly when $`\mu _i=ϵ_0`$. Instead, the resonance is shifted to lower $`|V|`$, by an amount $`T`$, so that it shifts with $`T`$. This is due to charge accumulation on the particle as the Fermi function sweeps by the energy of the quantum state, which limits the current at larger $`|V|`$ on account of Coulomb blockade. Our observations of this effect (not shown) are equivalent to results of Deshpande et al. , although the data in  were compared to an approximation that differed from Eq. (1) .

Non-equilibrium Regime: We can controllably tune the device so that more than one quantum state can participate in tunneling. This is illustrated by following line II in Fig. 1(b). This line corresponds to processes which are initiated by an electron tunneling off the nanoparticle from the quantum state 0 to electrode L. However, as one follows line II to higher V, past negative-sloping resonance lines which intersect line II, these lines indicate that the subsequent tunneling of an electron from electrode R back onto the nanoparticle can proceed via many different energy levels other than state 0. The total current under these conditions can be modeled by a master equation which takes into account all allowed transitions between the energetically-accessible $`n_0`$\- and ($`n_01`$)-electron states . In Fig. 2, we show tunneling diagrams depicting representative accessible states for resonances I-IV in a non-equilibrium regime.

In Fig. 3(a) we plot the step height in current associated with resonance lines I-IV, as $`V_g`$ and $`V`$ are tuned to follow the lines in the $`|V|`$-$`V_g`$ plane. Peaks I and III have approximately constant amplitude, while the currents for peaks II and IV grow quickly as $`|V|`$ enters the non-equilibrium regime. This can be understood trivially. For peaks I and III, the tunneling threshold is across the higher-resistance tunnel junction R. This junction is always rate-limiting and it matters little how many transport channels are available across junction L. For peaks II and IV, the tunneling threshold is across the low-resistance barrier L, but the rate-limiting process is how quickly electrons can tunnel across the other barrier. As $`|V|`$ is increased, more quantum levels contribute to this process, and the current grows. By measuring the current along peaks II and IV as levels are added one by one, and then fitting to master-equation results, we can measure rate-limiting tunneling rates for each quantum state: the $`\gamma _{iR}`$, for $`i=0\mathrm{to}5`$, are $`(5.3,15.7,8.0,16,15,9)\times 10^7s^1`$, $`\overline{\gamma _L}3\times 10^9s^1`$, and (assuming relaxation effects are negligible) $`\gamma _{1R}17\times 10^7s^1`$.

Figure 3(b) shows the widths of the tunneling resonance lines I-IV. These were determined by fitting the $`V`$ dependence of each conductance peak to the derivative of a Fermi function, and then converting the voltage width to effective temperature by multiplying by the capacitance ratio $`(e/k_B)C_L/(C_L+C_R)`$ for peaks I and III or $`(e/k_B)C_R/(C_L+C_R)`$ for peaks II and IV . In either the equilibrium or non-equilibrium regimes, the prediction of the simplest master equation  is that the peak shape should be a derivative of the Fermi function with a width approximately equal to the electronic temperature in the electrodes. Our measurements agree with this model within the equilibrium regime, with a constant electron temperature $`T`$ 90 mK. This is significantly higher than the $`T`$ achieved in non-gated tunneling devices using our apparatus, 45 mK, and we ascribe the difference to heating by leakage current from the gate to electrode R. Line III is broader than the others at low $`V`$ because the magnetic field of $`60`$ mT applied to drive the Al electrodes normal produces an unresolved Zeeman splitting ($`\mathrm{\Delta }E/k_B=2\mu _BH/k_B=80`$ mK).

As $`|V|`$ is increased into the non-equilibrium regime, peaks I and III undergo large increases in width and peak IV broadens slightly, while peak II shows no measurable change. The differences are not merely an effect of heating in the electrodes, because peaks II and IV have the largest magnitudes of current and power. We suggest that these measurements can be explained as a consequence of electronic interactions in the non-equilibrium regime, by a mechanism due to Agam et al. . Consider resonance line III, for which the tunneling threshold corresponds to an electron entering quantum state 0. For $`V>0.50`$ mV the next tunneling event, which discharges the particle, may occur out of different, lower energy states (see Fig. 2), leaving an electron-hole excitation on the nanoparticle. Agam et al. suggested that if this non-equilibrium state does not relax before the next electron tunnels onto the particle, it can shift the energy of the tunneling resonances on account of an alteration of the electron-electron interaction energy . In past work on smaller aluminum particles, shifted transitions were resolved individually ; however the relative shift is expected to decrease with increasing nanoparticle size , so it is reasonable that the shifts would produce only broadened resonances for the 10 nm particle under investigation here. Because a growing ensemble of different non-equilibrium states can be excited with increasing $`|V|`$, this mechanism can explain the increase in width of line III as a function of $`|V|`$. The same non-equilibrium mechanisms should also come into play for line IV, for $`V<0.64`$ mV, but the broadening here is reduced because the threshold tunneling event is across the lower-resistance junction, L. Barrier R quickly becomes rate limiting as line IV is crossed, so that higher-energy non-equilibrium resonances do not add significant additional current. In order for the non-equilibrium mechanism to apply for line III, the relaxation rate of some non-equilibrium excitations to the ground state must be comparable to or slower than $`\gamma _{0R}=5.3\times 10^7s^1`$. The rate predicted by Agam et al. for spin-preserving energy relaxation in aluminum particles is $`10^8s^1`$ .

Resonances I and II are a different case, because the tunneling threshold corresponds to an electron leaving quantum state 0. The subsequent tunneling event, adding an electron back to the nanoparticle, may for large $`V`$ occur into higher-energy states (see Fig. 2(a,b)), but nevertheless this excitation alone cannot produce a non-equilibrium shift in the energy of subsequent discharging transitions. The reason is that only this electron is free to tunnel off the nanoparticle; there is no electron in quantum state 0 whose transition energy might be shifted. Therefore within the picture of Agam et al. , no non-equilibrium broadening should be expected for levels I and II, in conflict with the data for level I. This discrepancy can be explained if non-equilibrium excitations on the nanoparticle can be generated not only by the tunneling transitions on or off the particle that have been considered previously , but also by transitions in which a high-energy electron relaxes within the nanoparticle and produces an electron-hole excitation. Broadening would then be generated by the Agam mechanism. Within this scenario, the difference between the broadening visible for resonance line I and the lack of broadening of line II would follow from the fact that for peak II a high-energy electron on the particle can quickly exit through the low-resistance tunnel junction L, while for peak I the high-energy particle must exit through the high-resistance junction R, giving a much longer residence time during which relaxation transitions can occur. In order for line I to be broadened, the fastest relaxation rates must become comparable to $`\gamma _{0R}=5.3\times 10^7s^1`$ as $`|V|`$ increases.

By tuning $`V_g`$ and $`V`$ into the non-equilibrium regime, the apparent energies of the $`dI/dV`$ peaks can also be changed. This is clearest for line III (Fig. 3(c)), which undergoes a shift of $`33\mu \mathrm{V}`$ to lower voltage when the threshold for non-equilibrium tunneling via state -1 (line V in Fig. 1(b)) is crossed. Because we have measured rate-limiting tunneling rates for the energetically-accessible states from the current amplitudes, we can test whether this shift can be explained by the simplest master equation , which assumes that the underlying energies of the quantum states are not changed by non-equilibrium interactions. Only one relevant parameter is not determined previously: $`x=\gamma _{0L}/\gamma _{1L}`$. The solution of the master equation does predict a voltage shift ($`T`$) for the conductance peak compared to the equilibrium case, and for $`x<0.15`$ it can explain the full value of the experimental shift. However, we judge this to be improbable, because the measured values of $`\gamma _{iR}`$ fall within a more narrow distribution. For $`x1`$ in Eq. (2), the predicted shift is $`15\mu \mathrm{V}`$ – much smaller than we measure. We can more naturally explain the full value of the shift by again taking into account that the presence of a non-equilibrium excitation can change the energy of a tunneling transition. In the equilibrium regime, the occupation of quantum state 0 corresponds to a transition from a fully-paired superconducting state on the aluminum particle to a state with one high-energy quasiparticle; in the non-equilibrium case, the transition can be from a state with two quasiparticles to one, with a transition energy lowered by $`2\mathrm{\Delta }`$ 0.35 meV . This is much bigger than $`k_BT10\mu \mathrm{eV}`$, and in this case the observed resonance is shifted to lower $`|V|`$ by an amount $`T`$ because electrons in the tail of the Fermi distribution can excite the non-equilibrium state and open the lower-energy current channel. For $`x=1`$ and $`T90`$ mK the master-equation result is that the measured shift can be produced by a non-equilibrium lowering of the transition energy by any amount greater than $`20\mu \mathrm{eV}2k_BT`$.

In summary, we have made a systematic study of the transition energies, widths, and current levels for electrons tunneling through discrete quantum states in a nanometer-scale electronic device. When transport occurs through a single quantum level, the results are in agreement with the expectations of sequential tunneling. At large values of $`|V|`$, the non-equilibrium population of excited electronic states, together with electron-electron interactions, can modify the widths and apparent energies of the tunneling resonances.

Acknowledgements: NSF DMR-0071631, the Packard Foundation, and the Cornell Nanofabrication Facility.


text:
# Static NLO susceptibilities: testing approximation schemes against exact results

## I Introduction

Photonic and optoelectronic applications heavily rely on materials with high nonlinear optical responses, and the development of new and more efficient materials for nonlinear optics (NLO) is a key issue in present days research. In this context, organic materials play an important role in view of their large and fast NLO responses and of the tunability of their properties via chemical synthesis. In order to trace guide-lines to organic synthesis, much effort is devoted to understand the origin of nonlinearity in these systems and to relate NLO responses to electronic structure and molecular geometry. Mobile electrons are required for large nonlinearity and conjugated molecules and/or polymers are materials of choice for NLO applications. Vibrations are strongly coupled to conjugated electrons: the effect of electron-phonon (e-ph) coupling in NLO responses at optical frequencies is not fully understood yet, but general agreement has emerged on the importance of e-ph coupling in static NLO responses.

Static susceptibilities are the opposite of the successive derivatives of the ground state (gs) energy with respect to a static electric field ($`F`$):

$$=_0\mu F\frac{1}{2}\alpha F^2\frac{1}{6}\beta F^3\frac{1}{24}\gamma F^4+\mathrm{}$$
(1)
where $`_0`$ is the gs energy at zero field, $`\mu `$ is the permanent dipole moment, $`\alpha `$ is the linear polarizability, and $`\beta `$ and $`\gamma `$ are the first and second hyperpolarizabilities, respectively. As it is always the case for energy derivatives, two different methods can be adopted to evaluate susceptibilities. The first method, the so-called finite-field (FF) approach, basically relies on direct numerical derivation techniques. The relevant system is the molecule in the presence of an external static electric field. The corresponding Hamiltonian is diagonalized to calculate the gs energy (or dipole moment) at different $`F`$ values, and (hyper)polarizabilities are calculated as numerical derivatives of the energy (or dipole moment) with respect to the electric field. The second approach treats $`F`$ as a perturbation on the molecular states, leading to the well-known sum over states (SOS) formulas. The main advantage of SOS over FF is that only unperturbed ($`F=0`$) eigenstates are required. The main disadvantage is that all eigenstates and transition dipole moments are required, whereas FF only requires the gs energy (or dipole moment). Of course the two approaches are perfectly equivalent if exact energies and eigenstates, including vibrational contributions, are inserted into equations. But for molecules and, more generally, for molecular materials, exact vibronic eigenstates are hardly accessible, and several approximations are introduced. Approximated calculations based on FF or SOS approaches are clearly different and it is important to investigate the merit and reliability of the various approximation schemes.

The first approximation usually introduced is the adiabatic (Born-Oppennheimer, BO) approximation. This is an almost ubiquitous approximation in molecular and condensed matter physics: accounting for the different masses of electrons and nuclei (and then for the different time scales of their motions), BO approximation separates electronic and vibrational degrees of freedom, so that each state is factorized into a product of an electronic and a vibrational wavefunction. This approximation is fairly safe: deviations from the adiabatic regime are expected only for degenerate or quasi-degenerate states, for narrow-band solids, or near to phase transitions.

In SOS approaches BO factorization is used to separate electronic and vibrational contributions to susceptibilities. Electronic contributions are those where all summations run on vibronic states in the excited electronic manifolds, whereas vibrational contributions have at least one of the indices in the summation running over excited vibrational states within the gs manifold. In the same spirit as BO, clamped nucleus (CN) approximation has been proposed to calculate SOS expressions for susceptibilities in a two-step procedure. In the first step the electronic (hyper)polarizabilities are calculated with nuclei clamped at their equilibrium positions. The resulting electronic contributions to (hyper)polarizabilities only contain sums over purely electronic states. In the second step the nuclear motion on the gs potential energy surface (PES) is accounted for, and vibrational contributions to susceptibilities are obtained. The adiabatic approximation, as directly applied to SOS or within the CN formalism, makes the calculation of susceptibilities possible, provided that electronic and vibrational eigenstates are known. However, in view of the anharmonicity of the relevant PES, calculating the vibrational states in all electronic manifolds (as required in BO-SOS formalism) or just in the gs manifold (as required in CN-SOS) is not easy. Often SOS calculations of static (hyper)polarizabilities of molecules and/or polymers are carried out in the harmonic approximation, where the solution of the vibrational problem can be written in closed form. Harmonic schemes work in the hypothesis that the electronic system responds linearly to the nuclear motion: harmonic approximation is bound to fail in systems where vibrations are coupled to mobile electrons, i.e. to electrons characterized by large nonlinear responses. In fact large anharmonic corrections to NLO responses have been recognized in recent ab initio calculations on conjugated molecules and polymers.

When applied to FF approaches, BO quite naturally leads to a two-step procedure. First of all the adiabatic electronic Hamiltonian (i.e. the total Hamiltonian minus the nuclear kinetic energy, KE, term) is solved to get the gs electronic energy, that represents the potential energy (PE) for the nuclear motion. Its derivatives on $`F`$ measure the PE contributions to susceptibilities. The calculation of the PE contribution is easy since only the gs electronic energy is required. Moreover, if nuclei are allowed to relax in response to $`F`$, susceptibilities calculated as PE derivatives naturally include vibrational contributions, and fully account for anharmonicity. The second step accounts for KE contribution: this is a more difficult step because the energy of the lowest vibrational state has to be calculated as a function of $`F`$, a non-trivial task due to anharmonicity. Moreover, in systems with large e-ph coupling, like conjugated polymers and/or oligomers, it can be difficult to get reliable calculations of $`F`$-dependent geometries, and numerical instabilities of the energy derivatives are sometimes observed.

In this paper we critically review the different approximation schemes for the calculation of static NLO responses, when applied to molecules of interest for NLO applications. In these molecules, in fact, electron-phonon coupling has particularly large and non-trivial effects, shedding doubts on the reliability of approximations that work well in more conventional systems. In particular we will focus on an interesting toy-model for mobile electrons coupled to vibrations: the Holstein $`DA`$ dimer. The so-called $`DA`$ dimer is a two-state model originally proposed by Mulliken to describe electron donor-acceptor charge transfer (CT) complexes in solution. Later on it was applied by Oudar and Chemla to describe push-pull chromophores where the donor and acceptor moieties are linked by a $`\pi `$-conjugated bridge to form a molecule. The extension of the model to include Holstein coupling to vibrations was first suggested by Person et al. to describe vibrational spectra of CT complexes, and was extensively discussed in the 80’s to describe many properties (ranging from vibrational spectra to structural instabilities) of CT salts. Holstein coupling in the $`DA`$ dimer has also been adopted to describe vibrations and/or polar solvation effects in push-pull chromophores: an extensive comparison with spectral data for several key molecules proved that this model contains the main ingredients to understand nonlinearity in these systems. In spite of its richness, the model is very simple and in the next Section we will describe its exact, non-adiabatic eigenstates, as well as its BO solution. Sections III and IV discuss the approximations usually introduced in SOS and FF approaches, respectively; Section V summarizes main results.

## II The Holstein $`DA`$ dimer: non-adiabatic and adiabatic eigenstates

The $`DA`$ dimer describes two electronic states, $`|DA`$ and $`|D^+A^{}`$, linearly coupled to a harmonic phonon, as depicted in Fig. 1. The relevant Hamiltonian is ($`\mathrm{}=1`$):

$$=2z_0\widehat{\rho }\sqrt{2}t\widehat{\sigma }_x+\frac{1}{2}\left(\omega ^2Q^2+P^2\right)\sqrt{2ϵ_{sp}}\omega Q\widehat{\rho }\widehat{\mu }F.$$
(2)
The first two terms in the above equation describe the electronic Hamiltonian, where $`\widehat{\rho }=(1\widehat{\sigma }_z)/2`$ is the ionicity operator, and $`\widehat{\sigma }_x`$, $`\widehat{\sigma }_z`$ represent the Pauli spin operators. The third and fourth terms describe a harmonic vibration of frequency $`\omega `$ ($`Q`$ and $`P`$ are the corresponding coordinate and momentum), coupled to the electronic system. The strength of the coupling is measured by $`ϵ_{sp}`$, the relaxation energy of $`|D^+A^{}`$ (see Fig. 1a). The last term accounts for a static electric field, $`F`$, interacting with the dipole moment operator, defined, by $`\widehat{\mu }=\mu _0\widehat{\rho }`$, with $`\mu _0=D^+A^{}|\widehat{\mu }|D^+A^{}`$. In the following, energies are expressed in $`\sqrt{2}t`$ units. We underline that the single-mode Hamiltonian in the above equation captures the essential physics of e-ph coupling, at least if one is not interested in the detailed description of vibrational spectra. The extension to the multimode case, possible in principle, is very demanding if the non-adiabatic solution is required.

The above Hamiltonian can be numerically diagonalized on the basis of the direct product of the two electronic states $`|DA`$ and $`|D^+A^{}`$ and of the reference vibrational states (i.e. the eigenstates of the harmonic oscillator in the third term of Eq. (2)). The basis is truncated by fixing a maximum number of phonon states, $`M`$; the corresponding $`2M\times 2M`$ matrix can be diagonalized up to fairly large $`M`$ values, yielding numerically exact non-adiabatic eigenstates. The minimum $`M`$ required to get convergence depends on the model parameters and on the properties of interest. All results presented in the following have been obtained with $`M=50`$.

The exact non-adiabatic eigenstates at $`F=0`$ are inserted into standard SOS expressions for static susceptibilities to get the exact NLO responses reported as continuous lines in Fig. 2 for $`ϵ_{sp}=1`$ and a few $`\omega `$ values. By the way, exactly the same results are obtained, within FF approach, from the successive derivatives of the exact non-adiabatic gs energy vs the applied field. In the same figure, dot-dashed lines correspond to the bare electronic susceptibilities, $`\chi _0^{(n)}`$, i.e. to the response of the two-state model with no e-ph coupling. The deviations of the continuous lines from the dot-dashed lines measure the vibrational contribution to static responses. Dashed lines show the PE contribution to susceptibilities. Analytical expressions for these curves have already been reported in Ref. : PE-susceptibilities are independent of $`\omega `$. Fig. 2 clearly shows the evolution of static susceptibilities with phonon frequency: in the low-$`\omega `$ limit, KE contributions vanish, and the exact curves tend to the limiting PE results. Vibrational contributions to static susceptibilities are very large in this limit, and increase with the order of nonlinearity. With increasing $`\omega `$, the vibrational contributions to static NLO responses decrease: in the antiadiabatic limit ($`\omega \mathrm{}`$) the vibrational contributions to static hyperpolarizabilities vanishes, and the exact curves tend to the bare electronic responses.

For real molecules or complexes, the large number of electronic states makes non-adiabatic calculations very demanding, and BO approximation is usually invoked. In BO approximation the effective electronic Hamiltonian, $`_{el}=1/2P^2`$, is defined by subtracting the nuclear KE term from the total Hamiltonian. In our case $`_{el}`$ describes two electronic states separated by an energy gap that linearly depends on $`Q`$, according to $`2z(Q)=2z_0\sqrt{2ϵ_{sp}}\omega Q`$. By diagonalizing the electronic $`2\times 2`$ matrix one gets analytical expressions for the $`Q`$-dependent ground and excited state energies, as shown in Fig. 1b. It is important to recognize that, even if the Hamiltonian in Eq. (2) is defined in terms of harmonic reference states (Fig. 1a), as a consequence of e-ph coupling, the ground and excited state PES are anharmonic (Fig. 1b). The anharmonicity of the potential prevents the analytical solution of the vibrational problem on either the ground or excited state; however numerically exact vibrational states can be calculated in both manifolds. In particular, the eigenstates of the harmonic oscillator with frequency $`\omega `$, centered at the relevant equilibrium position, are a good basis for the vibrational problem on either PES. The corresponding vibrational Hamiltonian is the sum of a KE term, whose matrix elements are trivial in the adopted basis, plus a PE term, whose matrix elements are calculated via numerical integration. Of course the vibrational matrix is diagonalized on a basis truncated to a large enough number of phonon states as to get convergence. Results presented in this paper have typically been obtained with 20 phonon states. Once BO eigenstates are obtained, the (transition) dipole moments entering SOS expressions can be calculated via numerical integration.

Static susceptibilities calculated within BO approximation are indistinguishable (in the scale of Fig. 2) from the exact ones, as long as $`\omega 0.2`$. Fig. 3 compares non-adiabatic (continuous lines) and BO (dashed lines) static susceptibilities for $`\omega =0.5`$, where deviations appear. Of course BO approximation becomes worst with increasing $`\omega `$ and is totally untenable for $`\omega 1`$. The adopted value of the small polaron binding energy, $`ϵ_{sp}\sqrt{2}t`$ applies to both push-pull chromophores and to CT complexes and salts. In push-pull chromophores the typical value of $`\sqrt{2}t1`$ eV is much larger than typical vibrational frequencies ($`\omega 1000`$ cm<sup>-1</sup>) and the dimensionless $`\omega 0.10.2`$ is safely in the BO regime. In CT salts, instead, $`\sqrt{2}t`$ is much smaller ($`0.2`$ eV) and the corresponding dimensionless $`\omega 0.51`$ suggests seizable non-adiabatic effects in NLO responses of these (narrow-band) systems.

## III SOS susceptibilities

The separation of electronic and vibrational degrees of freedom in BO approximation not only leads to simpler calculations, but also allows for additional insight. When applied to SOS expressions, BO quite immediately leads to a separation of electronic and vibrational contributions. If $`|Rv`$ indicates a BO vibronic state, i.e. the product between an electronic state $`|R`$ and a vibrational state $`|v`$ in the $`R`$ manifold, the SOS expressions for static susceptibilities read:

$`\alpha `$ $`=`$ $`2{\displaystyle \underset{R,v}{}}{\displaystyle \frac{0G|\widehat{\mu }|RvvR|\widehat{\mu }|G0}{\omega _{Rv}}}`$ (3)
$`\beta `$ $`=`$ $`6{\displaystyle \underset{R,v,S,u}{}}{\displaystyle \frac{0G|\widehat{\mu }|RvvR|\overline{\mu }|SuuS|\widehat{\mu }|G0}{\omega _{Rv}\omega _{Su}}}`$ (4)
$`\gamma `$ $`=`$ (6)
$`24`$ $`{\displaystyle \underset{R,v,S,u,T,w}{}}{\displaystyle \frac{0G|\widehat{\mu }|RvvR|\overline{\mu }|SuuS|\overline{\mu }|TwwT|\widehat{\mu }|G0}{\omega _{Rv}\omega _{Su}\omega _{Tw}}}`$ (7)
$``$ $`24{\displaystyle \underset{S,u,T,w}{}}{\displaystyle \frac{0G|\widehat{\mu }|SuuS|\widehat{\mu }|G00G|\widehat{\mu }|TwwT|\widehat{\mu }|G0}{\omega _{Su}^2\omega _{Tw}}},`$ (8)
where the gs, $`|G0`$, is excluded from the summations, $`\omega _{Rv}`$ is the frequency of the $`|G0|Rv`$ transition, and $`\overline{\mu }=\widehat{\mu }0G|\widehat{\mu }|G0`$. In the two-state model, the summations on $`R`$, $`S`$ and $`T`$ only run on the two electronic states $`|G`$ and $`|E`$. Then, according to the standard definition, the electronic contribution to susceptibilities is given by terms in Eq. (4) with $`R=S=T=E`$, and the vibrational contribution is described by terms where at least one of $`R`$, $`S`$ or $`T`$ states coincides with $`|G`$. We underline that electronic contributions to BO-SOS susceptibilities do not coincide with the responses of the bare electronic system (dot-dashed lines in Fig. 2). In fact summations in Eq. (4) run on vibronic states and not on pure electronic states. In the reliability range of BO approximation, closure on vibrational states in the excited state manifold works fine in reducing sums on intermediate vibronic states into sums on pure electronic states, but the property is evaluated for the true vibronic gs ($`|G0`$) and not as a pure electronic property (mediated over $`|G`$): the electronic contribution to BO-SOS susceptibilities partly include zero-point vibrational average (ZPVA) corrections, as defined in Ref..

The CN approximation introduces a different way of partitioning electronic and vibrational degrees of freedom. In this approach the electronic contribution to the $`n`$-th order susceptibility is calculated as the $`n`$-th order bare electron susceptibility, $`\chi _0^{(n)}`$. Two more contributions are then added to get the total response: (a) ZPVA corrections, that account for the difference between the electronic susceptibility calculated at the bottom of the gs PES (i.e. at the equilibrium geometry) and that relevant to the ground vibronic state (due to anharmonicity the relevant geometry does not coincide with the equilibrium one). These corrections are usually calculated via an expansion of electronic susceptibilities on $`Q`$; (b) vibrational contributions, that are calculated in terms of summations running on the vibrational states in the gs manifold only, as detailed in Ref..

BO-SOS and CN-SOS are based on the same adiabatic approximation and lead to basically identical results at least as long as the adiabatic approximation itself is reliable (approximately $`\omega <0.5`$). However, the partitioning of the susceptibilities is different and there is not a one-to-one correspondence between terms appearing in the two schemes: CN-SOS kills a few terms appearing in BO-SOS, and, at the same time, some of the terms, that are usually considered as corresponding, do not have exactly the same meaning in the two approaches. As a matter of fact, differences between CN and BO results amount to $`7\%`$ and $`15\%`$ in the electronic $`\beta `$ and $`\gamma `$, respectively, and $`2\%`$ and $`3\%`$ in the vibrational contributions to $`\beta `$ and $`\gamma `$ (these estimates refer to the parameters $`ϵ_{sp}=1`$, $`\omega =0.2`$; for $`\omega =0.1`$ deviations are even smaller). These small differences compensate each other, leading to negligible ($``$ 1%) differences in the total (electronic+vibrational) responses.

Our numerical solution of the adiabatic problem for the two-state model allowed us to test BO and CN approximations without introducing additional approximations. However getting exact vibrational eigenstates is hardly possible for complex (molecular) structures. For this reason CN-SOS calculations are often implemented by invoking the harmonic approximation for the vibrational problem. As already stressed in the Introduction, the harmonic approximation fails in systems where vibrations are coupled to electrons with large nonlinear responses. The relation between anharmonicity and nonlinearity is very well apparent in the adopted model, where both $`Q`$ and $`F`$ variables are coupled to the same electronic operator, $`\widehat{\rho }`$ (cf Eq. (2)). Successive $`Q`$-derivatives of the ground and excited state PE ($`_G`$ and $`_E`$, respectively) are therefore directly related to $`F`$-derivatives, i.e. to the bare electronic susceptibilities, $`\chi _0^{(n)}`$, as follows:

$$\frac{^n_{G/E}}{Q^n}=\omega ^2\delta _{2n}(2ϵ_{sp})^{n/2}\omega ^n\frac{\chi _0^{(n1)}}{\mu _0^n},$$
(9)
where $`\delta _{2n}`$ is the Kronecker-$`\delta `$, equal to 1 if $`n=2`$, zero otherwise. From this equation it turns out that systems with large hyperpolarizabilities, i.e. with large $`\chi _0^{(n)}`$, for $`n2`$, are characterized by largely anharmonic ground and excited state PES. As it has been discussed in Ref., the anharmonicity of vibrations coupled to delocalized electrons is hardly detected in coherent spectral measurements (e.g. electronic and vibrational absorption and/or Raman spectra), where experimental data can be satisfactorily reproduced in terms of parabolic PES with effective curvatures, but it shows up with large effects in incoherent spectral measurements (e.g. steady-state and/or time-resolved emission). Large anharmonic corrections are also expected in static nonlinear optical responses, as it is confirmed by extensive ab initio results on conjugated materials.

The dotted lines in Fig. 2 show the ($`\omega `$-independent) susceptibilities calculated in the best harmonic approximation (BHA). In this approach the ground and excited state PES are approximated by the two parabolas that best fit the two anharmonic PES at the gs equilibrium position (i.e. at the minimum of the gs PES): the relevant curvatures are then defined by Eq. (9) in terms of the bare electronic linear polarizability evaluated at the equilibrium geometry. Moreover, all electronic properties appearing in CN-SOS expressions are truncated to the linear term in $`Q`$. We defer a detailed comparison of BHA and exact curves to the next section, here we only underline that, as far as $`\alpha `$ is concerned, BHA and PE estimate exactly coincide, confirming that in the $`\omega 0`$ regime anharmonic corrections to the linear polarizability tend to vanish. In the same limit, nonlinear susceptibilities ($`\beta `$ and $`\gamma `$) are largely amplified by the anharmonicity of the gs PES, as demonstrated in panels (a) and (b) of Fig. 2, by the large deviations of BHA results (dotted lines) from PE curves (dashed lines).

## IV Finite-field susceptibilities: potential and kinetic energy contributions

In SOS approaches, susceptibilities are naturally partitioned into electronic and vibrational contributions. Within FF approach, on the opposite, PE contributions to the susceptibilities are naturally separated from contributions due to the nuclear KE. As already discussed in the Introduction, the lowest eigenstate of the electronic Hamiltonian as defined in the adiabatic approximation, is the PE for the motion of nuclei (in the gs manifold, of course). The total gs energy is obtained by summing the nuclear KE to the PE. Then, as long as BO applies, susceptibilities, i.e. the successive derivatives of the gs energy with respect to an applied electric field, can be calculated as sums of PE and KE $`F`$-derivatives. As a matter of fact, one can further separate the PE contribution to susceptibilities into an electronic and a vibrational part by simply comparing derivatives calculated at fixed nuclear geometry with those taken by allowing nuclei to relax following the application of the electric field. This separation is however fairly artificial and does not add so much to our understanding.

The nuclear KE vanishes in the $`\omega =0`$ limit, and the ($`\omega `$-independent) PE susceptibilities, reported as dashed lines in Fig. 2, represent the zero-frequency limit of the exact susceptibilities. Analytical expressions for the PE-susceptibilities of the Holstein DA dimer were already reported in Ref. . More generally, the calculation of PE-susceptibilities is easily implemented in quantum chemistry calculations, in fact it only requires the gs energy calculated at the relaxed geometry for different values of an externally applied field. Of course, geometry optimization is a crucial step in FF approach to hyperpolarizabilities and it is extremely important to fully relax the geometry for the chosen molecule in order to avoid spurious ‘strain’ or ‘instability’ contributions to hyperpolarizabilities.

The calculation of KE contributions is more difficult, since the $`F`$-dependence of the lowest vibrational state in the anharmonic gs PES is needed. In the adopted model, however, the calculation is feasible and, as long as BO applies, leads to basically exact results. Specifically, much as it occurs for BO-SOS and CN-SOS, for $`\omega 0.2`$ exact (continuous) curves in Fig. 2 are indistinguishable from BO-FF results. The nuclear KE contributes to susceptibilities in two different ways. First of all, due to anharmonicity, the molecular geometry in the vibronic gs is different from the equilibrium geometry (corresponding to the minimum of the gs PES). The PE susceptibilities have, in principle, to be corrected to account for this effect. The correction is however very small, and we found it always negligible in the investigated parameter range. The second contribution stems from the $`F`$-dependence of the nuclear KE itself: it is this contribution that indeed accounts for the deviations of the exact curves form the (dashed) PE curve (at least in the BO regime, $`\omega 0.2`$, where non-adiabatic corrections are negligible). KE contributions are of course very small for low $`\omega `$, but they increase with increasing $`\omega `$, leading to an overall decrease of the vibrational amplification of the static NLO responses. This is by no means accidental: in the antiadiabatic limit ($`\omega \sqrt{2}t`$) phonons cannot contribute to static susceptibilities and, with increasing $`\omega `$, KE contributions progressively increase to counterbalance PE contributions.

KE contributions to susceptibilities exactly vanish in the harmonic approximation. In fact, the equilibrium position in any harmonic vibrational state coincides with the bottom of the PES; moreover the nuclear KE is proportional to the harmonic frequency, i.e. to the curvature of the PES. For a parabolic PES, this quantity is obviously independent of $`Q`$, and hence of $`F`$ (cf Eq. (9)). We are then in the position to prove that the harmonic approximation is totally uncontrolled for the calculation of hyperpolarizabilities in any frequency regime. At low $`\omega `$, where PE contributions dominate, the harmonic approximation cannot account for anharmonic PES: this has no effect on the linear response, and the BHA estimate of $`\alpha `$ coincides with the exact curve in the $`\omega =0`$ limit. Anharmonic contributions to PE-hyperpolarizabilities are instead large, as demonstrated by the large deviations of BHA curves from the $`\omega =0`$ limiting curves. At the same time, not accounting for KE corrections, harmonic results do not represent a reliable approximation in the large $`\omega `$ regime either. It is possible that in a given parameter range BHA gives good estimates of a specific response, but this is due to an accidental compensation of two large errors. The harmonic approximation is unreliable in the calculation of hyperpolarizabilities, and, as it turns out clearly from Fig. 2, it is not possible to devise any parameter range where harmonic results represent a good approximation for all static susceptibilities.

## V Conclusions

The Mulliken $`DA`$ dimer model, extended to account for the Holstein coupling to vibrations (or equivalently to other slow degrees of freedom, including the orientational solvation coordinate), contains the essential physics to understand the subtle interplay between electrons, vibrations and external perturbations (namely an electric field) governing NLO responses of molecular materials. This simple toy-model is then very useful to test the reliability of several approximation schemes usually introduced in more refined quantum chemical description of these materials. It is important to recognize that some very well-known and very widely applied approximations can fail when applied to model the properties of materials with strongly nonlinear behavior. Perturbation theory, and particularly linear perturbation theory, is clearly inadequate in this context, and we already proved its inadequacy in reproducing spectral properties of push-pull chromophores as well as their static susceptibilities. More subtle effects have been demonstrated in two-photon absorption (TPA) spectra, where the standard Condon approximation badly fails even when it properly reproduces one-photon absorption spectra. Even more interestingly, an important vibrational channel is found to contribute to TPA spectra, in addition to standard electronic channels, leading to large effects in observed spectra. The harmonic approximations works fine in TPA spectra and more generally in vertical (coherent) processes, but it is inadequate to describe incoherent processes.

In this paper we have addressed the approximations usually introduced in the calculation of static NLO responses. The adiabatic approximation is fairly safe, at least as long as phonon frequencies are not too high (approximately $`\omega <0.5t`$, cf Fig. 3). It can be applied in different ways, BO or CN, in SOS approaches leading to a slightly different partition of electronic and vibrational contributions to susceptibilities, but to the same overall result. In FF it leads to a partitioning of potential and kinetic energy contributions to susceptibilities, but otherwise leads to the same result as in SOS. Closure over vibrational states is also a very good approximation, as long as the adiabatic approximation applies.

Potential energy contributions are clearly dominant over kinetic energy contributions in the low-vibrational frequency regime. The anharmonicity of the gs PES therefore gives very large corrections to hyperpolarizabilities: the harmonic approximation is untenable in this regime. With increasing phonon frequency, kinetic energy contributions globally reduce vibrational corrections until, in the antiadiabatic limit ($`\omega \mathrm{}`$) only electrons respond to static fields to regain purely electronic susceptibilities. Kinetic energy corrections vanish in the harmonic approximation and this approximation fails in the high frequency regime too.

Both FF and SOS approaches can be easily implemented in quantum chemistry calculations: the former has the advantage of only requiring knowledge about the lowest eigenstate, whereas all eigenstates are required in SOS. However only $`F=0`$ eigenstates are needed in SOS, whereas the $`F`$-dependence of the gs energy (or dipole moment) is required in FF. Anharmonicity is important and must be accounted for: in the low-vibrational frequency limit, where kinetic energy contributions are negligible, FF calculation is very convenient and fully accounts for anharmonicity. In the intermediate frequency regime the PE-FF calculation overestimates the vibrational amplification of NLO responses: either KE corrections have to be introduced in the FF scheme, or one has to make resort to SOS calculations, but in any case anharmonicity has to be accounted in order to get reliable estimates of hyperpolarizabilities.

ACKNOWLEDGEMENT

Work supported by the Italian National Research Council (CNR) within its “Progetto Finalizzato Materiali Speciali per Tecnologie Avanzate II” and by the Ministry of University of Scientific and Technological Research (MURST).


text:
# 1 Introduction

## 1 Introduction

A fruitful generalization of the duality between $`𝒩=4`$ $`SU(N)`$ gauge theory and type IIB strings on $`AdS_5\times S^5`$  is to consider other backgrounds of type IIB string theory, say on $`AdS_5\times X_5`$, where $`X_5`$ is a positively curved Einstein manifold. We are interested in theories with reduced supersymmetry, variety of gauge groups and matter contents.

The simplest example realizing the $`𝒩=1`$ super symmetry is provided by D3 branes at the singularity of a Calabi-Yau threefold with singularity known as conifold . The conifold in the AdS/CFT context was first considered in  and generalized to the case where conifold is deformed by the quantum effect of the fractional branes by Klebanov and Strassler .

To construct more realistic model, gauge theory with $`SO(N)/Sp(N)`$ is necessary and this can be obtained in brane language by including the orientifolds . Since there are so many $`Z_2`$ symmetry that can act on the conifold , one needs care to determine which symmetry is relevant to the specific orbifolding operation to define the orientifold. One way to fix the notion of the orientifold operation is to try to derive it from that of the type IIA brane construction where things are canonically defined. Under the T-duality the regular D4 branes are mapped to the regular D3 branes and and the fractional D4 branes are mapped to fractional branes which can be considerred as D5 branes wrapping the the vanishing a 2-cycle. We use $`\mathrm{T}`$ duality using the prescription given in .

There are discussions on orientifolding the conifold based on O6 plane in type IIA picture . Here we present a discussion on the orientifold in the conifold based on O4 brane of type IIA, because that is the one relevant to the physics of fractional branes discussed in Ref.. However it should be kept in mind that the resulting O3 branes has nature of fractional brane, since in type IIA picture O4 branes are between NS branes like D4 branes that is mapped to fractinal D3 branes. Therefore our resulting O3 branes has the character of the wrapped O5 branes, although our starting point in type IIA is completely different from the O6 brane configurations.

We will find that both $`O3^+`$ and $`O3^{}`$ are at the tip of the cone without annihilating each other since they are ’topologically protected’ from mutual annihilation. But their effect is combine to give zero net untwisted RR charges. The non-perturbative quantum effect deforms the tip of the cone. We will discuss how orientifold response to the deformation or resolution of the conifold. The finitely resolved conifold does not allow the orientifold, while deformed conifold leaves us an unresolved issue on supersymmetry. We will also show that we need two ‘fractional’ branes for the conformal invariance so that the theory should have gauge group $`SO(N+2)\times Sp(N)`$. We will give several explanations for this. We then give generalizations of the gauge theory results of ref.  to the case with orientifold. We describe the duality cascade and chiral symmetry breaking as well as the quantum deformation of the moduli space of the field theory in the presence of the conifold. We argue that the gravity solution is the same as Klebanov and Tsyetlin since it cannot distinguish the orientifolds and D branes.

The rest of paper comes in following order. In section 2, we give a brief review on the relevant background. In section 3, we derive the orientifold operation from that of Type IIA picture. In section 4, we perform renormalization group analysis to fix the that the conformally invariant configuration as well as to compare bulk and boundary theory. We also discuss the duality cascade and chiral symmetry breaking in the presence of the conifold. In section 5, we show that quantum moduli space of the corresponding gauge theory is a deformed conifold. We conclude in section 6.

## 2 Conifold and its type IIA brane construction: a review

A conifold is a complex submanifold in $`𝐂^\mathrm{𝟒}`$ described by the quadratic equation: $`_{i=1}^4z_i^2=0.`$ Its metric is known  to be $`ds^2=dr^2+r^2ds_{T^{1,1}}^2`$ with

$`ds_{T^{1,1}}^2={\displaystyle \frac{1}{9}}(d\psi +\mathrm{cos}\theta _1d\varphi _1+\mathrm{cos}\theta _2d\varphi _2)^2+{\displaystyle \frac{1}{6}}{\displaystyle \underset{i=1}{\overset{2}{}}}(d\theta _i^2+\mathrm{sin}^2\theta _id\varphi _i^2).`$ (1)
The base of the conifold, $`T^{1,1}`$, is an $`S^1`$ bundle over $`S^2\times S^2`$, and has the metric $`ds_{T^{1,1}}^2`$ given above. The conifold is a Calabi-Yau manifold and respects $`𝒩=2`$ supersymmetry. By putting $`N`$ D3 branes at the tip of the conifold we can get $`𝒩=1`$ supersymmetric field theory living on the D3. This field theory was constructed by Klebanov and Witten . It is a $`SU(N)\times SU(N)`$ gauge theory 111The reason for the gauge group being a product group is due to the underlying orbifold symmetry $`Z_2`$. coupled to two chiral superfields $`A_i,i=1,2`$ in the $`(N,\overline{N})`$ and $`B_j`$ in the $`(\overline{N},N)`$ representation. The most general superpotential which preserves the symmetry of the conifold is

$`W=\mathrm{Tr}ϵ^{ij}ϵ^{kl}A_iB_kA_jB_l.`$ (2)
One can wrap various D branes over the cycles of $`T^{1,1}`$ and identify these with states in the field theory .

A very intuitive way to understand the field theory content is to T-dualize the theory and consider brane configurations in type IIA string theory . As we discussed before, there are ambiguity which T-duality should we take. We focus the regular D3 branes and take $`\mathrm{T}_6`$ duality. In Ref. , $`N`$ D3 branes on a conifold is shown to be $`T_6`$-dual to the type IIA brane configuration with NS5(1,2,3,4,5), NS$`5^{}`$(1,2,3,7,8,9) and $`N`$ D4(1,2,3,6). $`x^6=\psi `$ is periodic, NS brane is at $`\psi =0`$, NS is at $`\psi =2\pi `$, $`N`$ D4 branes wrap the $`x^6`$ circle so that the system is an elliptic model. In type IIA picture it is very easy to see the resulting theory is $`SU(N)\times SU(N)`$ gauge theory with $`2N`$ flavors $`(A_i,B_i),i=1,2`$. The theory is conformally invariant, as can be checked explicitly from the calculation of beta functions .

Introducing $`M`$ fractional D3 branes into the picture, the resulting gauge theory is $`SU(N+M)\times SU(N)`$ (with gauge coupings $`g_1`$ and $`g_2`$ respectively) . In type IIA picture, this corresponds to putting $`M`$ D4 branes between $`0\psi 2\pi `$ and the theory is no longer conformally invariant. The two gauge couplings are determined as follows :

$$\frac{1}{g_1^2}+\frac{1}{g_2^2}e^\varphi ,\frac{1}{g_1^2}\frac{1}{g_2^2}e^\varphi \left[\left(_{S^2}B_2\right)1/2\right].$$
(3)
Since RG flow of the coupling constant in $`𝒩=1`$ gauge theory is logarithmic, the gravity dual is expected to have similar logarithmic behavior in the radial coordinate of $`AdS_5`$ space and in fact this is true . In terms of the type IIA brane construction, the two gauge couplings are determined by the positions of the NS5 branes along the $`x^6`$ circle. If one of the NS5 branes is located at $`x_6=0`$ and the other at $`x_6=a`$, then

$$\frac{1}{g_1^2}=\frac{l_6a}{g_s},\frac{1}{g_2^2}=\frac{a}{g_s},$$
(4)
where $`l_6`$ is the circumference of the $`x^6`$ circle .

As the NS5 branes approach each other, one of the couplings becomes strong. In fact, the two gauge couplings $`1/g_1^2`$ and $`1/g_2^2`$ flow in the opposite directions and there is a scale where one of the couplings diverge, which necessitates the use of Seiberg’s dual gauge theory . In the present situation, this corresponds to the moving the NS brane across the NS brane. After all the re-connections are made, we get $`SU(NM)\times SU(N)`$ theory . Notice that $`SU(N_fN_c)=SU(NM)`$.

As we go to the further IR region, the same process repeat until we get $`SU(M+p)\times SU(p)`$ with $`p`$ less than $`M`$. This is so called the cascade of the duality. The simplest case is $`p=1`$. one may consider this as one D3 brane probing the the background of M-fractional branes. Using the ADS superpotential , together with the original classical superpotetial for the conifold, it was shown that the conifold singularity is resolved into ‘deformed conifold’ . A very interesting phenomenon happens as a consequence, namely the chiral symmetry is broken.

## 3 Orientifold in conifold

The orientifolding operation is consist of world sheet orientation reversal and reflection of transverse spacetime, together with appropriate projection of the left moving fermion number. We will concentrate on the spacetime reflection part hereon. In the conifold there are many $`Z_2`$ operations  so that it is not clear which is the relevant $`Z_2`$ for the orientifolding. Therefore we want to derive the orientifold operation in the conifold from the type IIA picture where O4-plane is clearly defined.

We have D4 (01236), O4 (01236), NS5 (012345), NS5’(012389) with $`x^6`$ compact. Under the T-duality along the $`x^6`$, the whole configuration is mapped to the conifold. The prescription for the (T-duality) mapping flat $`R^6`$ to conifold suggested by Dasgupta and Mukhi  is:

$`x^4,x^5\mathrm{plane}`$ $``$ $`S^2\mathrm{described}\mathrm{by}\theta _1,\varphi _1,0\theta _i<\pi ,0\varphi _i<2\pi ,`$
$`x^8,x^9\mathrm{plane}`$ $``$ $`S^2\mathrm{described}\mathrm{by}\theta _2,\varphi _2,`$
$`x^6`$ $``$ $`\psi ,|x^7|\mathrm{log}1/r.`$ (5)
The orientifolding in type IIA picture is given by

$$x^ix^i\mathrm{for}i=4,5,7,8,9,\mathrm{and}x^ix^i,\mathrm{otherwise}.$$
(6)
Under T-duality along $`x^6`$, O4 brane becomes O3 brane. The reflections in $`x^4,x^5,x^6,x^8,x^9`$ induces an antipodal mapping in $`S^2`$ and $`S^3`$ of $`T^{1,1}`$, the base of the cone. Above reflection in terms of the polar coordinates is

$$\varphi _i\varphi _i+\pi _i,\theta _i\pi \theta _i,\psi \psi .$$
(7)
Now we want to express the above operation in terms of the conifold variables $`z_i`$’s. For doing that we have to express them in terms of the angular variables. Fortunately this has been done in Ref.. Introducing the variable

$$Z\left(\begin{array}{cc}z_{11}& z_{12}\\ z_{21}& z_{22}\end{array}\right)\frac{1}{\sqrt{2}}\left(\begin{array}{cc}z_3+iz_4& z_1iz_2\\ z_1+iz_2& z_3+iz_4\end{array}\right),$$
(8)
the equation defining the conifold can be rewritten as $`det(z_{ij})=0`$. The base $`T^{1,1}`$ is the intersection of the conifold with the $`S^5`$ sphere $`_{i=1}^4|z_i|^2=r^2`$, and is locally $`SU(2)\times SU(2)/U(1)`$. We parametrize $`T^{1,1}`$ by $`SU(2)`$ parameters $`a_i,b_i,i=1,2`$ :

$$\left(\begin{array}{c}a_i\\ b_i\end{array}\right)=\left(\begin{array}{c}\mathrm{cos}\frac{\theta _i}{2}e^{i(\psi +\varphi _i)/2}\\ \mathrm{sin}\frac{\theta _i}{2}e^{i(\psi \varphi _i)/2}\end{array}\right).$$
(9)
We can naturally map the $`S^2`$’s, which are compactifications of the 4,5 and 8,9 planes, into two $`SU(2)`$’s by the Hopf fibration $`\stackrel{}{n}=𝐚^{}\stackrel{}{\sigma }𝐚`$, where $`\stackrel{}{n}S^2`$, $`𝐚=\left(\begin{array}{c}a\\ b\end{array}\right)`$. With this identification, the antipodal mappings in the two $`S^2`$ induces a mapping in $`T^{1,1}`$ hence a mapping in the conifold. From the two $`SU(2)`$ matrices $`L=\left(\begin{array}{cc}a_1& \overline{b}_1\\ b_1& \overline{a}_1\end{array}\right),R=\left(\begin{array}{cc}a_2& \overline{b}_2\\ b_2& \overline{a}_2\end{array}\right)`$, we may construct the base variable

$$Z/r=LZ_0R^{}=\left(\begin{array}{cc}a_1b_2& a_1a_2\\ b_1b_2& b_1a_2\end{array}\right),$$
(10)
with $`r^2=\mathrm{Tr}(ZZ^{})=_i\left|z_i\right|^2`$.

In terms of $`SU(2)`$ variables the orientifold operation is

$$a_ii\overline{b}_i,b_ii\overline{a}_i,$$
(11)
which implies $`z_{11}\overline{z}_{22},z_{12}\overline{z}_{21}`$. In terms of original variables $`z_i`$’s, above operation has a simple expression:

$$z_i\overline{z}_i,i=1,\mathrm{},4.$$
(12)
The conifold is invariant under this $`Z_2`$ operation and so is the superpotential.

Let $`z_j=x_j+iy_j`$, $`j=1,\mathrm{},4`$ with $`x_j,y_j`$ be real. Then fixed points of this reflection are given by $`y_i=0,i=1,\mathrm{},4`$. Together with the conifold equation $`_iz_i^2=0`$, we conclude that $`z_i=0`$, i.e. the tip of the cone, is the only fixed point. This is consistent with the fact that the O3 is a point in the conifold.

However, O3 branes are like fractional D3 branes, since they are between NS branes rather than wrapping the whole circle of the IIA picture. In type IIB conifold picture, they are O5 branes wrapping the different singular $`S^2`$ cycles. In fact the base of the $`T^{1,1}`$ is $`\frac{SU(2)\times SU(2)}{U(1)}`$ and the $`U(1)`$ is acting symmetrically on both $`SU(2)`$’s  by

$$\left(\begin{array}{cc}e^{i\theta }& 0\\ 0& e^{i\theta }\end{array}\right).$$

It is also known that the second homology basis is a combination of two $`S^2`$’s above, namely, $`\mathrm{\Sigma }_2=S_2^2S_1^2`$ . The charge of the O3 brane is determined according to which $`S^2`$ the O5 wraps. Therefore both $`O3^+`$ and $`O3^{}`$ can be considered as a wrapped $`O5`$ branes wrapping different vanishing cycles. Since they are wrapping different $`S^2`$’s (although in the vanishing limit), we can say that they are stable due to the topological reason. However, since they are wrapping vanishing cycles at the same ’point’, they are like overlapping charges of opposite charge. So if we measure the net effect, there are no untwisted RR charges and it is an effect of one object that has only twisted RR charge. This situation is closely related to the orientifold of $`C^2/Z_N`$ considered by Uranga

We now ask what happen to the small resolution of the conifold?

$$\left|a_1\right|^2+\left|b_1\right|^2\left|a_2\right|^2\left|b_2\right|^2=\delta .$$
(13)
One immediately see that the small resolution of the conifold is invariant under the orientifold operation. However, from Eq.(11) there is no fixed point on the $`S^2`$ unless its size is zero. Only when the fractional branes are wrapping the ’vanishing cycle’, we can have orientifold in conifold. This means that the resolved conifold does not admit an super symmetric orientifold.

If there are large number of fractional branes, due to the quantum effect, the conifold background is modified  to deformed conifold $`_iz_i^2=\mu `$. In this case, the fixed points form a manifold $`S^3`$ given by $`_ix_i^2=\mu `$. Does this mean that O6 is created? Then, where is the orientifold in the deformed conifold? Part of the answer also lies in the large $`N`$ geometric transition : when the number of the fractional branes are large, it goes to the large N dual description where branes disappear and only flux remains. Presumably O3 branes disappeared leaving only its flux. Since O3- comes with extra two D3, $`O3^\pm `$ contribute the same amount of 3-form flux. Therefore the equality $`O3^+=O3^{}+2D3`$ holds as far as supergravity solutions are concerned. The flux of the fractional D3 brane charge resolve the singularity in the $`O3^+`$ side. Since O6 RR charge can not be created by the deformation process, the most natural answer to above question seems to be that the $`O3^+`$ charge is smeared uniformly over the fixed manifold, while $`O3^{}`$ is wrapping the vanishing $`S^2`$. We will be back to this issue in the discussion section for other possibilities. Now we make some remarks.

1. In type IIA, $`x^7`$ is a spectator. However, under the identification $`|x^7|=\mathrm{log}1/r`$, we are abandoning the region $`r1`$. We could equally cut out the $`r1`$ by the identification $`|x^7|=\mathrm{log}r`$. We have chosen above convention since we are interested in the near the singularity: $`r0`$ corresponds to $`x^7\mathrm{}`$.
2. Recently the same operation was considered in Ref. in the context of the topological string with O3 brane wrapping $`S^3`$. Here our interpretation is different: The fixed three sphere of the deformed conifold is due to smearing of the O3 brane charge.
3. It may be worth while to describe the orientifolding induced by the $`T_\varphi `$ duality, although we do not follow it. In this scheme, 4,5 and 8,9 plane is consequence of the T-duality along the circle action at the two sphere of the resolved conifold: two fixed points of the circle action under the $`T_\varphi `$ duality is mapped to the NS and NS branes. The orientifold operation of the type IIA branes does not involve the reflection along $`x^6=\psi `$ therefore not a antipodal point mapping. The action is simply given by

$$\varphi _i\varphi _i+\pi ,\theta _i\theta _i,\psi \psi .$$
(14)
In terms of $`a_i`$ and $`b_i`$, this is written

$$a_ia_i,b_iib_i.$$
(15)
Finally, in terms of the $`z_i`$’s,

$$z_1z_1,z_2z_2,z_3z_3,z_4z_4.$$
(16)
which also appear in literature as an orientifold operation.

## 4 The RG analysis and Duality cascade in the presence of the orientifold

If $`N`$ D3 branes are sitting at the singular point of the conifold, the gauge theory content turned out to be $`SU(N)\times SU(N)`$ . In this section, we discuss how introducing the orientifold modify the theory.

We use the type IIA brane picture. If we add O4 plane along 01236 direction, we have to change the sign of the RR charge of the O4 plane as we cross the NS brane . The corresponding gauge theory must be of alternating SO and Sp gauge groups. So one may naively expect that the gauge group would be $`SO(N)\times Sp(N)`$. (Here we use the convention where $`Sp(2)=SU(2)`$.) However, this is not conformal as we can see from the RG flow:

$`{\displaystyle \frac{d}{d\mathrm{log}(\mathrm{\Lambda }/\mu )}}\left[{\displaystyle \frac{8\pi ^2}{g_{so}^2}}\right]`$ $``$ $`3(N2)2N(1\gamma ),`$
$`{\displaystyle \frac{d}{d\mathrm{log}(\mathrm{\Lambda }/\mu )}}\left[{\displaystyle \frac{8\pi ^2}{g_{sp}^2}}\right]`$ $``$ $`3(N/2+1)(N)(1\gamma ).`$ (17)
One cannot require two beta functions to vanish simultaneously. We can overcome the difficulty if we replace $`SO(N)`$ by $`SO(N+2)`$. In brane language, we should add an extra brane and its image over the negative charged orientifold for the stable configuration.

The presence of the extra two branes can be also understood from the brane dynamics : for stable configuration, RR charge density must be continuous across the NS5. Since O4± has RR charge $`\pm 1`$, we need extra 2 D4 over the O4-. Then the gauge group is $`SO(N+2)\times Sp(N)`$. One may interpret the extra 2 branes as the fractional branes in IIB picture . So in the presence of the orientifold, the brane configuration with no fractional branes is not stable fixed point.

If we place $`N`$ D3 branes and $`M+2`$ fractional D3 branes on the conifold with an O3 plane, we obtain an $`SO(N+M+2)\times Sp(N)`$ or $`Sp(N+M)\times SO(N+2)`$ gauge groups. Here, we discuss the first one since the other one is exactly parallel. The two gauge group factors have holomorphic scales $`\mathrm{\Lambda }_1`$ and $`\stackrel{~}{\mathrm{\Lambda }_1}`$. The matter consists of two chiral superfields $`A_1,A_2`$ in the $`(𝐍+𝐌+\mathrm{𝟐},\overline{𝐍})`$ representation and two fields $`B_1,B_2`$ in the $`(𝐍+𝐌+\mathrm{𝟐},𝐍)`$ representation. The invariants of the theory under the global symmetry $`SU(2)\times SU(2)\times U(1)`$ are

$$I_1\lambda _1^{3M+2}\frac{\stackrel{~}{\mathrm{\Lambda }_1}^{2b_{Sp}}}{\mathrm{\Lambda }_1^{b_{SO}}}[\mathrm{tr}(A_iB_jA_kB_{\mathrm{}}ϵ^{ik}ϵ^j\mathrm{})]^{2M},$$
(18)
$$R_1^{(1)}=\frac{\mathrm{tr}[A_iB_j]\mathrm{tr}[A_kB_{\mathrm{}}]ϵ^{ik}ϵ^j\mathrm{}}{\mathrm{tr}(A_iB_jA_kB_{\mathrm{}}ϵ^{ik}ϵ^j\mathrm{})};$$
(19)
$$J_1\lambda _1^{b_{SO}+2b_{Sp}}\mathrm{\Lambda }_1^{b_{SO}}\stackrel{~}{\mathrm{\Lambda }_1}^{2b_{Sp}},$$
(20)
where $`b_{SO}`$, $`b_{Sp}`$ are $`\beta `$ function coefficients:

$`b_{SO}`$ $`=`$ $`3(N+M+22)4{\displaystyle \frac{N}{2}}1=N+3M,`$
$`b_{Sp}`$ $`=`$ $`3({\displaystyle \frac{N}{2}}+1)4{\displaystyle \frac{N+M+2}{2}}{\displaystyle \frac{1}{2}}=N/2M+1.`$ (21)
These statement can easily be verified by the quantum number assignment similar to that given in table 1 in Ref.. The only difference is that $`\mathrm{\Lambda }^{3(N+M)2N}`$ and $`\stackrel{~}{\mathrm{\Lambda }}^{3(N)2(N+M)}`$ of the for $`SU(N)\times SU(N)`$ should be replaced by $`\mathrm{\Lambda }^{3(N+M+22)2N}`$ and $`\stackrel{~}{\mathrm{\Lambda }}^{3(N+2)2(N+M+2)}`$ of the $`SO(N+M+2)\times Sp(N)`$. The superpotential of the model will get multiplicative renormalization depending on $`I_1,J_1,R_1`$’s.

In the presence of the orientifold, the geometry of the base of the cone is $`RP^2\times S^3`$. The 3-form flux breaks the conformal symmetry and $`B_2`$ gets radial dependence.

$$_{\mathrm{𝐑𝐏}^2}B_2(M/2)\mathrm{ln}(r/r_0),$$
(22)
which indicate the logarithmic running of $`\frac{1}{g_{SO}^2}\frac{1}{g_{Sp}^2}`$ in the $`SO(N+M+2)\times Sp(N)`$ gauge theory. To check this bulk result, we look at the $`\beta `$ functions of the boundary theory:

$`{\displaystyle \frac{d}{d\mathrm{log}(\mathrm{\Lambda }/\mu )}}\left[{\displaystyle \frac{8\pi ^2}{g_{so}^2}}\right]`$ $``$ $`3(N+M)2N(1\gamma )`$
$`{\displaystyle \frac{d}{d\mathrm{log}(\mathrm{\Lambda }/\mu )}}\left[{\displaystyle \frac{8\pi ^2}{g_{sp}^2}}\right]`$ $``$ $`3(N/2+1)(N+M+2)(1\gamma )`$ (23)
For the conformal invariance of $`M=0`$ case, we require $`\gamma =1/2`$. Notice that two flows give the same condition of the anomalous dimension in spite of the difference of the gauge group. The difference of the flow is

$$\frac{8\pi ^2}{g_{SO}^2}\frac{8\pi ^2}{g_{Sp}^2}=(4\gamma )M\mathrm{log}(\mathrm{\Lambda }/\mu ).$$
(24)
showing the the consistency of bulk and boundary result.

Since $`1/g_{SO}^2`$ and $`1/g_{Sp}^2`$ flow in opposite directions, there is a scale where the $`SO(N+M+2)`$ coupling, $`g_{SO}`$, diverges. Using the Seiberg duality transformation, we get the $`SO(2N[N+M+2]+4)=SO(NM+2)`$ gauge group with $`2N`$ flavors $`(a_i,b_i)`$ and “meson” bilinears $`M_{ij}=A_iB_j`$. The fields $`a_i`$ and $`b_i`$ are anti-fundamentals and fundamentals of $`Sp(N)`$, while the mesons are in the adjoint-plus-singlet of $`Sp(N)`$. The superpotential after the transformation

$$W=\lambda _1\mathrm{tr}M_{ij}M_k\mathrm{}ϵ^{ik}ϵ^j\mathrm{}F_1(I_1,J_1,R_1^{(s)})+\frac{1}{\mu }\mathrm{tr}M_{ij}a_ib_j,$$
(25)
where $`\mu `$ is the matching scale for the duality transformation , shows the $`M_{ij}`$ are actually massive. Thus we may integrate them out and get the superpotential

$$W=\lambda _2\mathrm{tr}a_ib_ja_kb_{\mathrm{}}ϵ^{ik}ϵ^j\mathrm{}F_2(I_2,J_2,R_2^{(s)}).$$
(26)
Here $`F_2`$, $`\lambda _2`$, $`I_2`$, $`J_2`$ and $`R_2`$ are defined similarly as in the original theory. Thus we obtain an $`Sp(N)\times SO(NM+2)`$ theory which resembles closely the theory we started with. This can be shown more carefully using the matching condition as discussed in . The next step is that the $`Sp(N)`$ gauge group becomes strongly coupled, and under a Seiberg duality transformation the full gauge group becomes $`SO(N3M+2)\times Sp(N4M)`$, and so forth.

## 5 Deformation of the orientifolded conifold

The classical field theory reveals that it represents branes moving on a conifold  by having the conifold as its moduli space of the gauge theory. Klebanov and Strassler showed that the vacuum configurations of gauge theory, the conifold, is modified to a deformed conifold by studying non-perturbative quantum corrections . Here we study corresponding phenomena in the presence of the orientifold. In our case, the minimal configuration is $`SO(M+4)\times Sp(2)`$, corresponding to a D-brane and its image in the presence of $`M+2`$ fractional branes. The fields are

$`A_{i,\alpha }^r,B_{j,\alpha }^a:i=1,2,\alpha =1,\mathrm{},N_c,r=1,\mathrm{},N_f.`$ (27)
Define $`N_{ij}`$ by $`N_{ij}^{rs}=_\alpha A_{i,\alpha }^rB_{j,\alpha }^s`$. Then the classical superpotential can be written as

$$W_c=\frac{1}{2}\lambda \text{Tr}A_iB_jA_kB_lϵ^{ik}ϵ^{jl}=\frac{1}{2}\lambda \text{Tr}N_{ij}N_{kl}ϵ^{ik}ϵ^{jl}.$$
(28)
where the trace is over the flavor index. For the gauge theory with $`N_c`$ colors and $`N_f`$ flavors, the quantum effect gives the ADS super potential . In our case, the gauge group is $`SO(N+M+2)\times Sp(N)`$. So one may expect that we should add up the ADS potential for each gauge group. However, one should notice that we are interested in $`M2=N`$ and the ADS potential is meaningful only for $`N_c>N_f`$. So we only have to consider the superpotential for the $`SO(N+M+2)`$ part. We propose that following is the leading order superpotential for the product group:

$$W_{total}=\lambda W_c+(N_c2N_f)\left(\frac{\mathrm{\Lambda }^{b_0}}{det_{ir;js}N_{ij}^{rs}}\right)^{1/(N_c2N_f)}.$$
(29)
Now we consider the completely higgsed configuration $`SU(N)U(1)^N`$, where the matrix $`N_{ij}`$ have vacuum expectation values:

$$N_{ij}^{rs}=n_{ij}^{(r)}\delta ^{rs}.$$
(30)
This is justified since $`r`$-th and $`s`$-th fractional (flavor) branes are far separated in higgsed case. So, the determinant is trivially factorized:

$$\text{det}_{ir,js}N_{ij}^{(rs)}=\underset{r=1}{\overset{N_f}{}}W^{(r)},$$
(31)
where $`W^{(r)}=n_{ij}^{(r)}n_{kl}^{(r)}ϵ^{ik}ϵ^{jl}`$ Using $`W_c=\lambda _rW^{(r)},`$ and considering $`N_f=2`$ case, the total superpotential becomes

$$W_{total}=\lambda (W^{(1)}+W^{(2)})+(N_c2N_f)\left(\frac{\mathrm{\Lambda }^{b_0}}{W^{(1)}W^{(2)}}\right)^{1/(N_c2N_f)},$$
(32)
where $`b_0=3(N_c2)N_f`$. The total potential is minimized if we have

$$(N_{ij}N_{kl}ϵ^{ik}ϵ^{jl})^{(N_cN_f)}=\left(\frac{\mathrm{\Lambda }^{b_0}}{\lambda ^{NcN_f2}}\right).$$
(33)
Notice that $`N_cN_f=(M+2)2=M`$ in our case. This is nothing but the equation for the deformed conifold and there are $`M`$ branches: each of the probe branes move on the deformed conifolds.

The M-theory curve  for the Type IIA version is given as follows:

$$(vw)^M=\mathrm{\Lambda }^{3M},t=v^{N_cN_f}=v^M.$$
(34)
The curve indicates that if there is more than two fractional branes ($`M>0`$), the whole brane configuration is nicely connected and this indicates the deformation of the conifold. Interesting fact is that the curve is insensitive to the presence of the regular D branes since it contribute to $`N_f`$ as well as to $`N_c`$. It depends only on the number of fractional branes. Orientifolds does not change the formal behavior at all, either. It just require that $`M`$ is even.

## 6 Conclusion

In this paper, we discussed the string theory in conifold in the presence of the orientifold. We first mapped the orientifolding operation from the type IIA brane picture to the conifold picture, using the T-duality along the $`x^6`$ direction. Under the T-duality the fractional D4 branes ( D4 branes between NS-NS branes) are mapped to fractional D3 branes, which is identified as D5 wrapping the vanishing 2-cycles. We showed how the conifold can admit orientifolds of both $`+`$ and $``$ charges. Blowing up the tip, the base is still $`S^2\times S^2\times S^1`$. The fractional $`D3`$ and $`O3^{}`$ are D5 and $`O5`$ branes respectively wrapping one of the $`S^2`$. The $`O3^+`$ is O5 brane wrapping the other $`S^2`$. Since homology 2-cycle which defines the RR charge is given by the difference of two $`S^2`$’s, it is natural to have both $`O3^+`$ and $`O3^{}`$ simultaneously at the tip of the cone without annihilating each other: they are topologically protected not to be annihilated. We then showed that in the presence of the orientifold, the conformally invariant configuration is $`SO(N+2)\times Sp(N)`$ rather than $`SO(N)\times Sp(N)`$. This is shown both by field theory and brane dynamics. If we add fractional branes, there are duality cascade as in Ref.. We analyzed the corresponding gauge theory as well as the super gravity solution. When there are fractional branes, the conifold is deformed as is the case of $`SU(N)\times SU(N)`$. We showed this by writing quantum corrected superpotential for the product group.

However more detailed discussion of the supergravity part is not discussed. We expect that the gravity solution is the same as Klebanov and Tsyetlin since it cannot distinguish the orientifolds and D branes. However we expect that there are minor modification due to torsion part. Also the K-theoretic discussion of the orientifold charge in the conifold is also worthwhile do be discussed in detail. We wish to come back to this issue in near future.

Finally we give a discussion on the deformation of the orientifold. As we have seen in section 6, the gauge theory analysis showed that the conifold is deformed by the quantum mechanical effect. First we should notice that, once the conifold is deformed, the connection between $`z_{ij}`$ and $`(a_i,b_i)`$ is lost. So there is no way to derive the orientifold operation in the deformed conifold from that of the IIA branes. What we have done is to assume that the orientifold operation of the original (singular) conifold $`z_i\overline{z}_i`$ is still valid for the deformed case. Then $`S^3`$ is the fixed points, which gave us a big conceptual trouble, since the most natural interpretation is to regard $`S^3`$ as the cycle wrapped by O6. But this is not allowed by the SUSY of IIB. In this situation, we have several options.

1. The identification of orientifold operation $`z_i\overline{z}_i`$, for the deformed case, is not correct.
2. O3 charge is smeared out so that it has O3 charges but have the effect of the O6. In this case, the SUSY is intact.
3. O6 really is created by the deformation. In this case we loose SUSY.

We excluded the first option since $`z_i\overline{z}_i`$ is the only one that leads to the correct one in the limit of zero deformation parameter. We have chosen the second option since gauge theory or SUGRA so far does not show any signal of supersymmetry breaking.

However, before the paper of Strassler and Klebanov  appeared, Mukhi and DasGupta in  claimed that SUSY is broken when the singularity is resolved. So the situation not entirely clear even in the absence of the orientifold and we believe that it is an interesting subject to study. One may further speculate that the appearance of fixed $`S^3`$ means appearance of O6 meaning the supersymmetry is broken. If true it can be used as a dynamical supersymmetry breaking mechanism. However, this require more extensive analysis which goes beyond the scope of present work and we hope to settle this issue in later publication.

Acknowledgements
We would like to thank I. Klebanov and Uranga for discussion and APCTP for support and hospitality during our visit. The work of SJS is supported by KOSEF 1999-2-112-003-5. The work of CA and SN is supported by KOSEF 2000-1-11200-001-3. This work is also supported by BK21 program of Korea Research Fund (2001).