text:
# A Measurement of the Flux of Cosmic Ray Iron at 5×1013 eV

## 1 Introduction

It is important to extend direct measurements of cosmic ray composition towards the knee of the cosmic ray spectrum. Determining composition in this region is one of the keys to understanding cosmic ray acceleration and propagation. Air shower techniques are a potential tool for measuring composition at the knee, but their interpretation ultimately relies on extrapolation of particle interaction models to regimes not directly tested in laboratory experiments. Emulsion techniques and electronic detectors must be improved before they can provide an independent calibration of air shower experiments, and so the potential for systematic error exists. We have designed an experiment (Evenson and Seckel 1995; Seckel et al. 1999) to determine the flux of iron nuclei at energies between 3.0$`\times 10^{13}`$ and 3.0$`\times 10^{14}`$ eV by detecting atmospheric Cherenkov radiation from the incoming nucleus before it undergoes a nuclear interaction. By achieving an absolute measurement of the iron flux, relatively free of the uncertainty in particle interaction models, the Balloon Air CHerenkov (BACH) technique can provide an independent calibration for air shower experiments, thus opening a path for composition measurements for energies at and above the knee. This technique was discussed by Sitte (1965) and Gough (1975), and used initially by Sood (1983; Sood and Panettieri, 1981). In this paper we report the result obtained from the initial flight of our developmental balloon payload (PDQ BACH).

## 2 Technique

A high energy iron nucleus entering the atmosphere has an interaction length of 13.2 gm/cm2 (Hufner 1985). Before interacting, a sufficiently high energy nucleus will begin to radiate photons via the Cherenkov process. These photons travel close to the primary particle, forming a “light pool” whose characteristics can be calculated from elementary considerations and knowledge of the density profile of the atmosphere. Figure 1 illustrates the radial profile of the light pool for a variety of particle energies and nuclei. At 35 km altitude, attainable by a NASA balloon, the diameter of the lightpool increases from zero at the threshold energy to a maximum of about 15 meters for a fully relativistic particle incident 45 degrees from the zenith. The intensity has a central peak, a cusp at the edge of the pool, and a broad floor at intermediate radii. For an iron nucleus the floor is about 2000 photons per square meter in the sprectral band from 0.3 to 0.6 microns. With finite size light collectors the peak and cusp are softened, and for the 24.3 cm radius light collectors employed by PDQ BACH the floor corresponds to about 400 photons. Quantum efficiency of our system is about 0.1, so a minimum of 40 photoelectrons are detected in each collector. These arrive with a time spread of less than 0.1 ns and produce a single pulse when observed by a photomultiplier tube. The trigger resolving time of the instrument is determined by the approximate 4 ns pulse width of the photomultipliers.

Cherenkov flashes from iron must be detected against the night sky background light. Fluctuations in response to background light can cause false event triggers reducing live time of the instrument, and sufficiently strong fluctuations may be confused with signals from cosmic ray iron nuclei. The PDQ BACH telemetry system takes 0.78 seconds to readout one event. To achieve livetime of 90% we must limit the background trigger rate to 0.1 Hz. This means that the false trigger probability should be less than $`4\times 10^{10}`$ per resolving time. Operating two collectors in coincidence allows a fluctuation trigger rate $`P=2\times 10^5`$ per resolving time in each collector.

The mean number of background photoelectrons in one sample is $`N_B=ϵ\varphi _BA\mathrm{\Omega }t`$, where $`ϵ`$ is the detection efficiency, $`\varphi _B`$ is the background flux, $`A`$ is the collector area, $`\mathrm{\Omega }`$ is the solid angle, and $`t`$ is the integration time. Pre-flight, we anticipated a dark sky background of $`10^8\mathrm{cm}^2\mathrm{sr}^1\mathrm{s}^1`$, which corresponds to $`N_B4`$ photoelectrons for BACH’s design parameters. For this $`N_B`$, a single collector threshold of 14 detected photoelectrons yields $`P=10^5`$. The equivalent photon count is shown in Figure 1 as the horizontal line labeled “P5”. For $`N_B=8`$ the corresponding threshold rises to 22 photoelectrons. During flight, $`N_B`$ varied between 6 and 8 photoelectrons, with occasional excursions as bright stars came into the field of view. At 40 photoelectrons, cosmic ray iron events lie well above the thresholds of 14-22 photoelectrons required from deadtime consideration. The fluctuations themselves are on the tail of a steep distribution, so an increase in threshold by a modest $`15\%`$ in post-flight analysis lowers the effective background rate to less than one per hour and allows for a clean separation of fluctuations and cosmic ray events.

Lower $`Z`$ nuclei may be confused with iron if the cusp or central peak of the light pool falls on a collector. By separating the two coincident detectors by $`3`$ m one may ensure that the central peak cannot simultaneously overlay the two collectors, thus largely eliminating the background from low $`Z`$ nuclei. The cusps for CNO nuclei are not strong enough to generate triggers, but silicon events where the cusp falls on both collectors may mimic iron. As an illustration of the characteristics of the BACH instrument, Figure 2 shows the results of exposing a numerical model of the instrument to a simulated flux of cosmic rays. Panel (a) shows that silicon may provide some contamination for low intensity iron events, while panel (b) demonstrates that CNO nuclei do not. Details of the model are given in Section V.

We also considered other sources of background, such as Cherenkov radiation from 10 GeV electrons entering a collector ($`50`$ photons) and 100 GeV iron nuclei generating flourescent radiation directly in front of the collectors ($`10`$ photons), but neither has enough intensity to trigger the instrument. Numerous lower energy charged particles may pass through a detector generating spurious signals. By running two collectors in coincidence, the trigger rate from such processes is diminished. Offline analysis can distinguish between the few remaining particle events and true cosmic ray events.

In evaluating the viabilty of the BACH technique one must also consider details of the atmospheric density profile and factors that could alter the assumed straight line trajectory of the radiating nucleus. We assume a US standard atmosphere (NOAA 1976) profile with the density normalized to match the density of the atmosphere at the pressure altitude of the balloon. We have not modeled absorption due to ozone as the edge of that feature is obscured by absorption in the optical components of the instrument. We considered deflection by multiple Coulomb scattering and hadronic interactions (both $`10^5`$ radians), as well as magnetic bending which may be a factor 5 larger depending on the track geometry. These effects are significantly smaller than $`10^3`$ radians, the Cherenkov angle at altitude for fully relativistic particles, and so have not been included in our simulations of the experiment. Hadronic interactions eventually disintegrate the primary nucleus; however, most of these interactions are soft. Disintegration takes several interaction lengths and occurs mostly after the primary has passed the balloon. We estimate a 9% exposure correction for this process, but do not model it in detail.

## 3 Instrumentation

Our payload is illustrated schematically in Figure 3. Given the intrinsic 100 psec duration of the Cherenkov light flash, the primary approach to skylight background is to reduce signal integration time. To this end, we employ Hamamatsu R4143 three inch diameter photomultipliers with a risetime of 2.4 nanoseconds at our operating point of 1200V. We also reduce background by limiting the acceptance aperture with light collectors (208.8 cm long, 24.3 cm radius aperture) based on concepts developed by Winston and colleagues (Welford and Winston 1989). These collectors have a nearly flat response within a well defined 7.5 degree half angle field of view. To reduce backgrounds, two collectors, separated on axis by 2.87 m, were operated in coincidence. The Winston collectors for PDQ BACH were constructed using a precision-machined aluminum mandrel. Each collector consisted of two upper and two lower sections, molded individually on the mandrel from carbon fiber epoxy composite material. A reflective aluminum coating was evaporated onto the sections, which were then assembled.

Photomultipliers attached to the collectors must operate in the presence of a comparatively high DC current from starlight, which must be monitored accurately while still maintaining satisfactory high frequency response. By operating with negative high voltage, the anode current can be directly measured; however, a photocathode (and, consequently, photomultiplier faceplate) maintained at negative high voltage while exposed to low atmospheric pressure poses an unacceptable risk of corona discharge. Our solution was to construct Pressurized Optical Detectors (PODs), as illustrated in Figure 4, to operate the photomultipliers at one atmosphere, and couple them to the collectors via light pipes of Bicron UV transparent acrylic. The light pipes extend out from the PODs to allow thermal isolation of the PODs and collectors, while providing electrical insulation as well. Surfaces of each light pipe are polished and kept from contacting the POD structure except for a 0.32 cm tapered flange that provides a pressure seal. Each light pipe is internally reflecting for almost all signal photon trajectories. The PODs also contain R329 photomultipliers coupled to plastic scintillators that cover the detector photomultipliers and internal sections of the light pipes to serve as an anticoincidence shield against cosmic rays. Each photomultiplier receives high voltage from an individual power supply, adjustable in flight.

Data were acquired with a Tektronics TDS640A digital oscilloscope which flew as part of the payload. This oscilliscope captures 2 gigasamples per second in four channels simultaneously. The four channels were used to record data from the two R4143 photomultipliers (denoted C1 and C2) and the two anticoincidence detectors (A1 and A2). Internal oscilloscope logic was used to detect coincident triggers where the signals in C1 and C2 both exceed thresholds. The thresholds were independently adjustable during flight. The anticoincidence channels were not part of the trigger, but were used in off line analysis. The TDS640A comes equipped with a GPIB interface, so we developed a microcontroller-based unit to interface it to the command and telemetry system.

Crucial to the scientific return of a BACH payload is calibration of photomultiplier gain, pulse shape, and the efficiency of the light collectors. The full optical assemblies were calibrated for angular efficiency and alignment within the Ft. Sumner hanger (Seckel, et al 1999). A 2.54 cm diameter diffuse light source was moved horizontally and vertically across the field of view at a distance of 7.5 m from the front of the collectors. The results of these tests were interpreted with the aid of ray tracing models for the collectors/PODs. The models include losses due to reflection from the collector (4%), absorption by the aluminum coating (20%), absorption in the light pipe (12%), losses from the optical path (4%), and yield an overall collector efficiency of 0.60 for a head on Cherenkov spectrum weighted by the photomultiplier quantum efficiency. Off axis, performance drops near the edge of the field of view in a calculable and testable manner. Minor adjustments were made to the individual collector models based on the test results. The adjusted models for the individual collectors were then used to determine instrument response to a light source at infinity.

Calibration of the photomultiplier response was complicated by low gain at the flight voltage of 1200 V and an inability to resolve pulses from single photoelectrons. We therefore followed a multistep process including determination of the absolute single photoelectron response at high voltage (2150 V), the relative gain curve from 2150 V down to 1200 V, and a comparison of charge and pulse height response at 1200 V. The mean pulse height per photoelectron for strong pulses at flight voltage was 0.24 mV/pe in collector C1 and 0.19 mV/pe in collector C2. Single photoelectron pulse height distributions were determined at high voltage and extrapolated down to flight voltage using an 8 stage model for the R4143 photomultipliers.

Ground level observations of cosmic rays provided a check on the absolute normalization of the collector and POD response. Muons passing through the light pipe serve as a rough check of the light pipe model and photomultiplier gain. Similarly, we are able to measure Cherenkov light from ground level muons and electrons passing through the collector volume on trajectories nearly parallel to the optical axis. We also performed ground level observations of Cherenkov light produced by cosmic ray air showers as a full system check of the payload in flight configuration.

## 4 Observations

PDQ BACH was launched at sunset from Ft. Sumner, NM, on April 22, 1998 and cut down $`6`$ hours later, yielding a useful exposure of 2.75 hours. Atmospheric pressure during the observation varied from 8.05 to 8.5 mb, corresponding to a threshold energy for iron of $`2.2\times 10^{13}`$ eV. The instrument collected 1793 triggers, each consisting of 250 ns (500 samples) of waveform from four channels: the two collectors and the two anticoincidence detectors. Since the signals resulting from the Cherenkov light are only a few nanoseconds wide, most of the waveform data are available for monitoring the skylight background. Expecting only a handful of true iron events, the trigger thresholds on the detectors were intentionally set to acquire events at a rate of approximately 0.1 Hz to allow detailed study of background sources.

Sample events from flight data are shown in Figure 5. Panel (a) contains one of the events identified as due to an iron nucleus. Note the identical arrival times of the similar sized pulses from the two collectors (C1 and C2) and the absence of signal from either anticoincidence (A1 and A2). Panel (b) shows a trigger with two moderate signals in C1 and C2. On the basis of pulse height, the event could be a silicon nucleus or iron incident at the edge of the field of view; however, there is a 2.5 ns offset in the two pulses, so this event is classified as a rare, high-amplitude skylight fluctuation. Panel (c) shows a typical event where a relativistic particle passes through one of the light pipes accompanied by a skylight fluctuation in the other channel. In addition to the signal in the anticoincidence detector, the event shows asymmetric response and a significant time lag between the C1 and C2 pulses. The time delay between collector (R4143) and anticoincidence (R329) signals results from the different transit times within the two types of photomultiplier. Panel (d) shows a rarer type of background event. The relative timing of the signals indicates that a particle passed through the light pipes and anticoincidence detectors in both of the PODs.

Data from one, 1.75 hour, segment of the flight are shown in Figure 6. This represents the longest exposure taken with unchanged trigger conditions. Data obtained under other sets of trigger conditions were analyzed separately, and the final results were combined. Several events from the 1.75 hour segment have the exact characteristics that we expect from cosmic ray iron nuclei, namely strong coincident sharp signals well above threshold in both channels and no signal in either anticoincidence detector. Of the remainder, the majority have amplitudes just above threshold, and are most likely due to simultaneous fluctuations in the skylight background. As expected, the rate of such triggers is a strong function of the amplitude of the DC current observed in the waveforms outside the pulse region. We observed numerous asymmetric events where C1 records a strong signal and C2 is just above threshold, resulting from charged particles passing through the light pipe in C1 accompanied by statistical fluctuations in C2. The rate of these light pipe events is also consistent with expectations. The relative lack of events where particle hits are evident in C2 coincident with a fluctuation in C1 results from an asymmetry in trigger levels on the oscilloscope. This was intentional, to aid in controlling the trigger rate in a changing background. About 75% of light pipe events are tagged by A1, consistent with our modeling for coverage by the anticoincidence detector. A few events have strong pulses in both channels but are clearly not caused by Cherenkov light, since both anticoincidence detectors also show signals. In all of these cases one signal leads the other by about 10 ns, consistent with a single charged particle passing through both light pipes.

## 5 Results

After eliminating events with anticoindence hits, the data can be described as the sum of two distributions characterized by the time separation $`dt`$ of the two pulses. For true cosmic ray events, the signals arrive within a narrow coincidence window corresponding to the geometry of the collectors ($`dt1`$ ns). Skylight fluctuations have a broader distribution corresponding to the pulse width of the photomultipliers ($`dt4`$ ns). Events with strong signals are predominantly of the first type, whereas events with smaller signals are mostly in the second. We count cosmic ray events as those with signal greater than 8.5 mV in both channels and $`dt<0.75`$ ns. There are 9 such events. Some of these events could be skylight events which just happen to have $`dt<0.75`$. By comparing to the clean skylight sample, i.e. events with $`dt>0.75`$ ns, we estimate a sklylight background count of 0.75 events. We also estimate that inefficiency of the anticoincidence detectors has resulted in false acceptance of 1 event, leaving a total of 7.25 cosmic ray events.

To convert this number to an iron flux, we compare to the mathematical model of the instrument exposed to a flux of cosmic ray nuclei of assumed composition and spectrum, as illustrated in Figure 2. The incident flux is taken to be isotropic over the field of view of the collectors with impact parameters that allow for all positions and orientations of the payload. Knowledge of the light pool, impact parameters, and POD/collector efficiency are combined to yield an expected number of photoelectrons, which in turn is the basis for drawing a signal amplitude from a distribution based on Poisson statistics, the photomultiplier single photoelectron response and the modeled pulseheight distributions for the photomultipliers. Contribution of background light to the signal is modeled by overlaying livetime weighted randomly chosen sections of waveform from the pre-trigger flight data. Pairs of pulseheights thus determined are displayed in the scatterplots of Figure 2. Applying the 8.5 mV cut to the simulation yields a count which is directly comparable to the 7.25 counts from analysis of the flight data.

Normalizing the exposure to the simulation, yields a flux of $`\varphi _{Fe}(5.05_{1.6}^{+4.5}\times 10^4\mathrm{GeV})=677_{139}^{+192}E^{2.5}`$ m-2 sr-1 s-1 GeV1.5. We use Poisson statistics to estimate the statistical error. The energy value (range) is the peak (width at half maximum) of the distribution obtained by convoluting the instrument response with an assumed $`E^{2.5}`$ spectrum. The result is reasonably independent of the 8.5 mV analysis threshold, but depends on the assumed iron spectral index. With a small number of events and uncertain calibration we cannot determine the spectrum. Our result also depends on the composition, as silicon nuclei with particular impact parameters can imitate the signals from iron. For the 8.5 mV threshold analysis, the result includes a 6$`\%`$ correction to account for mis-identified silicon. Corrections were also made for hadronic interactions of iron in the atmosphere (estimated 9 $`\%`$), and for contamination of the iron flux by nuclei with Z $`=`$ 17-25 (estimated 10$`\%`$).

Figure 7 compares our new determination of the iron flux with earlier measurements. Our results are consistent with both the CRN and JACEE determinations. All of the recent data are above the flux reported by Sood (1983). We have not reanalyzed Sood’s experiment, but have attempted to improve on his pioneering effort in several ways to make our instrument easier to model and the interpretation of the data more certain.

We thank A. McDermott, L. Piccirillo, E. Rode, C. Scharmberg, L. Shulman and J. Poirier for technical assistance. The collectors were produced in collaboration with the University of Bremen, the Astronomical Institute in Hamburg, and the private firm “Surf Affairs”, located near Kiel. The thank the National Scientific Balloon Facility for the balloon flight and the NSBF staff for preflight assistance. The work was supported in part by NASA grant NAG5-5063 and NAG5-5221.


format: text

text:
# 1 Introcuction

## 1 Introcuction

In year 2000, a number of improvements were made to the LEP collider. It increased Higgs discovery potential up to approximately 115 GeV if LEP combination were made. In LEPC held in November 2000, an excess of candidates, which is compatible with the Standard Model(SM) Higgs boson of 115 GeV, was reported . The results of a search for the SM Higgs boson with the OPAL detector at LEP was reported in this talk based on the results described in published letter .

## 2 The Analyses

The data used in the analyses correspond to integrated luminosities of approximately 216 pb-1 at 192–202 GeV, 80 pb-1 at 203–206 GeV, and 130 pb-1 at centre-of-mass energies ($`\sqrt{s}`$) higher than 206 GeV. During 2000 (1999), data were taken at $`\sqrt{s}=200`$-209 GeV (192-202 GeV) with a luminosity-weighted mean $`\sqrt{s}`$ of 206.1 (197.6) GeV.

The analyses follows our previous publication. The search was performed in the following final states, $`\mathrm{H}^0`$$`\mathrm{Z}^0`$$``$bb̄qq̄ (four-jet), $`\mathrm{H}^0`$$`\mathrm{Z}^0`$$``$bb̄$`\nu `$$`\overline{\nu }`$ (missing energy), $`\mathrm{H}^0`$$`\mathrm{Z}^0`$$``$bb̄$`\tau ^+`$$`\tau ^{}`$ and $`\tau ^+`$$`\tau ^{}`$qq̄ (tau), $`\mathrm{H}^0`$$`\mathrm{Z}^0`$$``$bb̄$`\mathrm{e}^+`$$`\mathrm{e}^{}`$ and $`\mathrm{H}^0`$$`\mathrm{Z}^0`$$``$bb̄$`\mu ^+`$$`\mu ^{}`$ (electron and muon).

In all channels, preselections followed by likelihood selections were applied to select the Higgs signal. The likelihood consists of b-tagging information, event shape variables, lepton identification, (in)compatibilities to the signal(background) mass hypothesis and so on. The likelihood selection were optimised using a mixture of various masses of signal Monte Carlo samples. Hence, the selection was not optimized to a particular mass.

The b-tagging variable $``$ is a likelihood combination of three independent b-tagging information, lifetime sensitive part, jet kinematics and lepton $`p_t`$ information. Lifetime sensitive part is an artificial neural network consists of four input variables including both displaced vertex and track impact parameter information. The performance of b-tagging was checked using the calibration data taken at $`\mathrm{Z}^0`$ peak in 2000, as well as higher energy data as shown in Figure 1.

The number of selected events in all search channels was 156 with 146.1$`\pm `$11.9 expected from the SM background processes. Figure 2 shows the likelihood distributions of each channels after preselection and the reconstructed mass distributions after likelihood selection.

After the event selections search channels were combined using the likelihood ratio method  or the method described in previous paper . Three different final discriminants, “test statistic”, were used in different channels to reflect the different background conditions. The four-jet channel uses a likelihood consists of four variables including b-tagging, reconstructed mass and two event shape parameters. The missing energy, electron and muon channel uses likelihood calculated from mass and selection likelihood output. The tau channel uses simply the mass distribution as a test statistic.

## 3 Results

Figure 3(a) shows the log-likelihood ratio as a function of test mass. The two bands in the figure are 68% and 95% contours centred on the median expectation. Figure 3(b) shows the projection of (a) at test mass of 115 GeV, corresponding to the probability density of $`2\mathrm{ln}𝒬`$ for signal+background and background only hypotheses. The background confidence level (1$``$CLb) and signal+background confidence level (CLs+b) are 0.2 and 0.4 shown in dark and light shaded area in the figure, respectively.

The Figure 4(a) shows 1$``$CLb as a function of the test mass. An excess can be seen in a range of 95–120 GeV with the lowest value of 0.02 at 107 GeV. The probability to observe such an excess anywhere in the range of test masses between 95 and 120 GeV is approximately 10%, estimated from the size of the range and the reconstructed mass resolution. The expected curve in the presence of 115 GeV Higgs is also seen in the figure.

The signal rate limit $`n_{95}`$ observed in data is shown as a function of $`m_{\mathrm{H}^0}`$ in Figure 4(b) together with its median expectation in background-only hypothesis. Figure 4 (b) also shows the expected accepted signal rate. A lower mass bound of 109.7 GeV was obtained, and the expected limit was 112.5 GeV.

## 4 Summary

A search for the Standard Model Higgs boson has been performed with the OPAL detector at LEP based on the full data sample collected at $`\sqrt{s}`$192–209 GeV in 1999 and 2000. The largest deviation with respect to the expected SM background in the background confidence level, 1$``$CLb, is 0.02 at 107 GeV. The observed excess is less significant than that expected for a SM Higgs boson with a 107 GeV mass. A lower bound of 109.7 GeV on the mass of the SM Higgs boson is obtained at the 95% confidence level while the median expectation for the background-only hypothesis is 112.5 GeV. For a Higgs boson with a mass of 115 GeV, 1$``$CLb is 0.2 while CLs+b is 0.4, indicating that the data slightly favour the hypothesis that a signal is present, but also that the data are consistent with the background hypothesis.


text:
# Dielectric Behavior of Nonspherical Cell Suspensions

## I Introduction

Dielectric spectroscopy has become a quantitative method of real-time monitoring of cell growth in suspensions \[1–3\]. The real-time monitoring has advantages over conventional techniques and would be important to investigate the dynamic behavior of cell growth. There are many factors that may influence the dielectric behavior of biological materials: structure, orientation of dipoles, surface conductances, membrane transport processes, etc. All these factors influence one another and it is difficult to separate out the effect of a single one. However, some effects can be dominant at certain ranges of frequencies. For instance, the dielectric dispersion spectrum of living cell suspensions in radiofrequencies is mainly determined by the cell shape. The objective of this work is to develop a theory for such correlation, on which new applications in biotechnology rely.

More recently, Asami  reported that the dielectric dispersion spectrum of fission yeast cells in a suspension was mainly composed of two sub-dispersions. The experimental data revealed that the low-frequency sub-dispersion depended on the cell length, while the high-frequency one was independent of it. Asami adopted a shell-ellipsoid model , in which an ellipsoid is covered with an insulating shell as the electrical model of nonspherical biological cells. The comparison between model calculation  and experimental data  was far from being satisfactory. Asami suggested that the discrepancy is attributed to the rod-like cell shape. However, the depolarization factor needed in his theory is difficult to calculate for cells of rod-like shape because of lack of available theories.

In this work, we propose the use of the spectral representation  for analyzing the cell models. The spectral representation is a rigorous mathematical formalism of the effective dielectric constant of a two-phase composite material . It offers the advantage of the separation of material parameters (namely the dielectric constant and conductivity) from the cell structure information, thus simplifying the study. From the spectral representation, one can readily derive the dielectric dispersion spectrum, with the dispersion strength as well as the characteristic frequency being explicitly expressed in terms of the structure parameters and the materials parameters of the cell suspension (see section III below). The actual shape of the real and imaginary parts of the permittivity over the relaxation region can be uniquely determined by the Debye relaxation spectrum, parametrized by the characteristic frequencies and the dispersion strengths. So, we can study the impact of these parameters on the dispersion spectrum directly.

## II Formalism

### A Spectral representation theory

We consider a two-component composite dielectric with complex dielectric constant $`ϵ(𝐫)`$ equal to $`ϵ_2=\epsilon _2+\sigma _2/j2\pi f`$ in the host medium and $`ϵ_1=\epsilon _1+\sigma _1/j2\pi f`$ in the embedding medium. An interface $`\mathrm{\Sigma }`$ separates the two media. In a uniform applied electric field $`𝐄=E_0\widehat{𝐳}`$ (for convenience, let $`E_0=1`$), the electrostatic potential satisfies the Laplace’s equation:

$$[(1\frac{1}{s}\eta (𝐫))\mathrm{\Phi }(𝐫)]=0,$$
(1)
where $`s=ϵ_2/(ϵ_2ϵ_1)`$ denotes the relevant material parameter and $`\eta (𝐫)`$ is the characteristic function of the composite, having value 1 for $`𝐫`$ in the embedding medium and 0 otherwise. The electric potential $`\mathrm{\Phi }(𝐫)`$ can be solved formally as:

$$\mathrm{\Phi }(𝐫)=\mathrm{\Phi }_0(𝐫)+\frac{1}{s}𝑑𝐫^{}\eta (𝐫^{})^{}G_0(𝐫𝐫^{})^{}\mathrm{\Phi }(𝐫^{})$$
(2)
where $`G_0(𝐫𝐫^{})=1/4\pi |𝐫𝐫^{}|`$ is the free space Green’s function, and $`\mathrm{\Phi }_0(𝐫)=𝐫\widehat{𝐳}=z`$ is the potential of the unperturbed uniform field E. It is instructive to convert the volume integration into the surface integration via the Green’s second identity and only deal with the potential on the interface $`\mathrm{\Sigma }`$ . We denote an integral-differential operator $`\mathrm{\Gamma }`$:

$$\mathrm{\Gamma }\mathrm{\Phi }(𝐫)=_\mathrm{\Sigma }^{}𝑑𝐒^{}^{}G_0(𝐫𝐫^{})\mathrm{\Phi }(𝐫^{})+\frac{1}{2}\mathrm{\Phi }(𝐫),𝐫\mathrm{\Sigma },$$
(3)
to avoid the singularity of $`G_0(𝐫𝐫^{})`$ when the integration variable $`𝐫^{}`$ approaches the point of $`𝐫`$ . The integration with a “prime” denotes the restricted integration which excludes $`𝐫^{}=𝐫`$. Let $`\mathrm{\Psi }_n(𝐫)`$ and $`s_n`$ be the $`n`$th eigenfunction and eigenvalue of the $`\mathrm{\Gamma }`$ operator respectively. We can expand $`\mathrm{\Phi }_0(𝐫)`$ and $`\mathrm{\Phi }(𝐫)`$ on the interface $`\mathrm{\Sigma }`$ in a series expansion of eigenfunction $`\mathrm{\Psi }_n(𝐫)`$:

$`\mathrm{\Phi }_0(𝐫)`$ $`=`$ $`{\displaystyle \underset{n}{}}z_n\mathrm{\Psi }_n(𝐫),`$ (4)
$`\mathrm{\Phi }(𝐫)`$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{sz_n}{ss_n}}\mathrm{\Psi }_n(𝐫),`$ (5)
where $`z_n`$ are the expansion coefficients. Then we can write the effective dielectric constant $`\overline{ϵ}`$ in the Bergman-Milton representation :

$`\overline{ϵ}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle 𝑑Vϵ(𝐫)E_z}`$ (6)
$`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle 𝑑Vϵ_2[1\frac{1}{s}\eta (𝐫)]\frac{\mathrm{\Phi }}{z}}`$ (7)
$`=`$ $`ϵ_2\left[1{\displaystyle \frac{1}{V}}{\displaystyle \underset{n}{}}{\displaystyle \frac{z_n}{ss_n}}{\displaystyle _\mathrm{\Sigma }}𝑑𝐒\widehat{𝐳}\mathrm{\Phi }_n(𝐫)\right]`$ (8)
$`=`$ $`ϵ_2\left[1p{\displaystyle \underset{n}{}}{\displaystyle \frac{f_n}{ss_n}}\right],`$ (9)
where $`p`$ is the volume fraction of the suspending cells. Note that $`E_z`$ is a dimensionless electric field because $`E_0=1`$. The eigenvalue $`s_n`$ and the spectral function $`f_n`$ can be proved to be real and satisfy simple properties that $`0<s_n<1`$, $`f_n>0`$ and $`f_n=1`$. We shall show that both the spectral function $`f_n`$ and the eigenvalue $`s_n`$ determine the dielectric behavior of cell suspensions.

### B Cells with an axis of revolution

Now the principal problem is to calculate the eigenvalue $`s_n`$ and the spectral function $`f_n`$. For many cells interacting with one another, it is a formidable task. However, in the limit of a dilute cell suspension and weak applied field, one can regard the cells in suspension as being noninteracting and randomly oriented and the problem is reduced to that of a single cell. We will consider cells with an axis of revolution, namely, the spheroidal and the spherocylinder cells  to mimic cells of rod-like shape. The prolate spheroid is generated by rotating an ellipse around its major axis, while the spherocylinder is obtained by fitting two hemispherical caps at both ends of a circular cylinder. For a prolate spheroid, the eigenvalues and eigenfunctions can be calculated exactly. The only nonzero $`f_n`$ equals unity for E being along the major or minor axis of the prolate spheroid, and the corresponding eigenvalues are given by:

$`s_z`$ $`=`$ $`{\displaystyle \frac{1}{q^21}}+{\displaystyle \frac{q\mathrm{ln}[q+(q^21)^{1/2}]}{(q^21)^{3/2}}},`$ (10)
$`s_x`$ $`=`$ $`(1s_z)/2.`$ (11)
where $`q`$ is the length to diameter ratio; $`z`$ and $`x`$ refer to the direction along the major and minor axis respectively. For a spherocylinder, the consideration of the symmetry properties of the cell will help us choose the appropriate orthogonal basis for calculating the matrix elements of the $`\mathrm{\Gamma }`$ operator. Because of the rotation symmetry about the major axis of the spherocylinder, the eigenfunction is necessarily of the form $`(a_n\mathrm{cos}n\theta +b_n\mathrm{sin}n\theta )f(x)`$ with $`n`$ being an integer. Due to the inversion symmetry of the cell, $`f(x)`$ must be either odd or even functions. It is convenient to expand $`f(x)`$ in a series of Legendre polynomials $`P_m(x/l)`$, where $`2l`$ is the length of the cell. The applied uniform field E can always be resolved into two components along the major and minor axes of the spherocylinder, so we can calculate the $`s_n`$ and $`f_n`$ for E along the major and minor axes separately. By symmetry, in order to obtain a nonzero $`f_n`$, the eigenfunction should be the form of $`A_mP_{2m+1}(x/l)`$ for E being along the major axis, while it reads $`\mathrm{cos}\theta B_mP_{2m}(x/l)`$ for E being along the minor axis, with $`m=0,1,2,\mathrm{}`$. Using this orthogonal basis, we can calculate a truncated matrix according to the precision needed. We should remark that the matrix is generally nonsymmetric.

## III Dielectric dispersion spectrum

We show here that from the spectral representation, one can readily derive the dielectric dispersion spectrum. Substituting $`ϵ_1=\epsilon _1+\sigma _1/j2\pi f`$ and $`ϵ_2=\epsilon _2+\sigma _2/j2\pi f`$ ($`\epsilon `$ and $`\sigma `$ being the real and imaginary parts of the complex dielectric constant) into Eq.(9), defining a new parameter $`t=\sigma _2/(\sigma _2\sigma _1)`$ and re-defining $`s=\epsilon _2/(\epsilon _2\epsilon _1)`$, we rewrite the effective dielectric constant $`\overline{ϵ}`$ after simple manipulations:

$$\overline{ϵ}=ϵ_H+\underset{n}{}\frac{\mathrm{\Delta }ϵ_n}{1+jf/f_n^c}+\frac{\sigma _L}{j2\pi f},$$
(12)
where $`ϵ_H`$ and $`\sigma _L`$ are the high-frequency dielectric constant and the low-frequency conductivity respectively, while $`\mathrm{\Delta }ϵ_n`$ are the dispersion magnitudes, $`f_n^c`$ are the characteristic frequencies of the $`n`$th sub-dispersion.

We have already shown that there are only two poles in the spectral representation of the prolate spheroids. In what follows, we will show that there are two dominant poles in the spectral representation of the spherocylinder and hence there are two sub-dispersions in the dielectric dispersion spectrum. The dispersion magnitudes and dispersion frequencies are given by:

$`\mathrm{\Delta }ϵ_1`$ $`=`$ $`{\displaystyle \frac{1}{3}}p\epsilon _2{\displaystyle \frac{s_1(st)^2}{s(ss_1)(ts_1)^2}},`$ (13)
$`\mathrm{\Delta }ϵ_2`$ $`=`$ $`{\displaystyle \frac{2}{3}}p\epsilon _2{\displaystyle \frac{s_2(ts)^2}{s(ss_2)(ts_2)^2}},`$ (14)
$`f_1^c`$ $`=`$ $`{\displaystyle \frac{\sigma _2s(ts_1)}{2\pi \epsilon _2t(ss_1)}},`$ (15)
$`f_2^c`$ $`=`$ $`{\displaystyle \frac{\sigma _2s(ts_2)}{2\pi \epsilon _2t(ss_2)}}.`$ (16)
Thus, we are able to obtain the dispersion strengths as well as the characteristic frequencies explicitly in terms of the structure parameters and the materials parameters of the cell suspension.

The dielectric dispersion spectrum of a dilute suspension of prolate spheroids is mainly composed of two sub-dispersions, namely, $`s_z`$ is responsible for the lower frequency one and $`s_x`$ for the higher one. For a spherocylinder, we obtain a nonvanishing series of $`f_n`$ and $`s_n`$. Along the major axis, $`f_1`$ is dominant for all $`q`$ and we can omit the smaller ones. This dominant $`f_1`$ is plotted in Fig.1 against $`q`$, and the corresponding $`s_z`$ are plotted in Fig.2, together with the exact result of a prolate spheroid. As is evident in Fig.1, we can see that the difference between the two models is indeed small. Along the minor axis, the solution becomes more complicated. The dominant $`f_2`$ near $`q=1`$ decreases quickly as $`q`$ increases; another $`f_3`$ increases and takes over at large $`q`$. These two $`f_n`$ are also plotted in Fig.1 and their corresponding eigenvalues are plotted in Fig.3. As shown in Fig.3, the two eigenvalues tend to that of a prolate spheroid in the limit of both small and large $`q`$.

Near $`q=2`$, the two $`f_n`$ have comparable values, resulting in two sub-dispersions at higher frequency. These sub-dispersions can interfere with each other, rendering it difficult to find the characteristic frequencies of the different sub-dispersions. Physically, the local field is the most nonuniform in this case. Nevertheless, we will consider cells of large length and omit this complication.

With Eqs.(13)–(16), it is easy to calculate the effect of the rod-like cell structure on the dispersion spectrum and to compare with experiment data. We will show that the spherocylinder model does give some improvement towards the experimental result. However, the improvement is too small to close up the gap as Asami expected. In fact, we shall see that Asami omitted the material parameters which will play an important role in the experimental condition. By introducing the conductivity contrast $`t=\sigma _2/(\sigma _2\sigma _1)`$, we found that a small negative $`t`$, i.e. $`\sigma _1\sigma _2`$, should be used to close up the discrepancy.

We estimate $`t`$ and $`s`$ by fitting Eqs.(13)–(16) to the experimental ratio of $`\mathrm{\Delta }ϵ_1/\mathrm{\Delta }ϵ_2`$ and $`f_2^c/f_1^c`$, and we get $`t=0.0014`$ and $`s=5.0`$. It means that $`\sigma _1700\sigma _2`$ and $`\epsilon _10.80\epsilon _2`$. The enhanced conductivity of cell cytoplasm is attributed to the membrane potential. The result is in contrast to the previous (unjustified) claim that $`\sigma _1\sigma _2`$.

Table I lists the $`\mathrm{\Delta }ϵ_1/\mathrm{\Delta }ϵ_2`$ ratio and $`f_2^c/f_1^c`$ ratio for both experimental and theoretical results. Using the fitting material parameters, the improvement is obvious for both the prolate spheroid model and the spherocylinder model, while the difference between the two models is quite small.

## Discussion and Conclusion

In this work, we have applied the spectral representation to the dielectric dispersion of suspensions of fission yeast cells. As mentioned by Asami , the discrepancies between theory and experiment may be attributed to the rod-like cell shape. For cells of nonconventional shape, however, there exists no available cell model in the literature and we must develop the spectral representation from first principles.

More precisely, we have developed a Green’s function formalism  for calculating the spectral representation of rods of finite length. We modelled the rod-like cells as the spherocylinders, i.e., circular cylinders with two hemispherical caps at both ends. We solved the spectral representation of the effective dielectric constant from first principles. Similar formalism was adopted for cell suspensions near their sub-division point \[9–11\].

Generally speaking, when the axial ratio $`q`$ is larger than 4, the prolate spheroid model can be employed as a good approximation for rod-like cell structures. For $`q<4`$, the dielectric behavior will become more sensitive to the cell structure of the suspending particles, and there are more sub-dispersions than that of the prolate spheroid suspensions.

Our model does not include the rotational or vibrational effects, and our results are expected to be valid only for weak electric fields. Our model also ignores the multi-shell nature of the cells. Usually the multi-shell model is used to explain the high-frequency steps of spherical cell suspensions. Similar conclusions were found in one of our previous paper on multi-shell dielectric spheres in electrorheological (ER) fluids  to account for the effects of water coating on the ER effects. In Ref., we also showed that the spectral representation can still be used for multi-shell model, albeit with a slight modification.

In Asami’s experiment, there exist three subdispersions, the highest frequency step (above 10MHz) is due to the vacuole and cell wall as mentioned by Asami, while the two lower frequency step is evidently dependent on the cell shape. And the disperson magnitude of the highest frequency step is much smaller than that of the two lower frequency ones. So one expects that the multi-shell model has only a small effect on the lower frequency steps. In fact the multi-shell model was used in Asami’s theory, but the discrepancy, as we mentioned in the text, is still large.

The large cytoplasmic conductivity is a key result of our investigation. We believe that the large cytoplasmic conductivity is reasonable because the cells have to maintain a higher ion concentration in their cytoplasm to avoid the shrinkage of cells due to a loss of water across the cell membrane. However, to our knowledge, there exists no direct experimental mesurement on the cytoplasmic conductivity. In our work, we propose a convenient and practical means of determining the cytoplasmic conductivity from the dielectric spectroscopy data, which analysis can be important for biotechnology.

## Acknowledgments

This work was supported in part by the Direct Grant for Research of the Research Committee and in part by the Research Grants Council of the Hong Kong SAR Government. K.W.Y. acknowledges useful conversation with Professor G. Q. Gu.

Fig.1/Lei, Wan, Yu, Sun

Fig.2/Lei, Wan, Yu, Sun

Fig.3/Lei, Wan, Yu, Sun


text:
# Nanoscopic stripe-like inhomogeneities and optical conductivity of doped cuprates

## 1 Introduction

High-$`T_c`$ superconductivity and other unconventional properties of doped quasi-2D cuprates remain a challenging and hot debated problem. As usual, the active $`CuO_2`$ planes in cuprates are considered as homogeneous systems, starting from one-, three- or multiband Hubbard models. Often only in-plane $`Cu(3d)O(2p\sigma )`$ basis set is taken into account and the well defined Zhang-Rice singlet ground state for doped holes is presumed. In frames of such theories the dependence of various properties of copper oxides on doping are related to the rearrangements of Hubbard bands.

However there are clear experimental evidences in favour of intrinsic and generic inhomogeneous nature of the systems under consideration with well developed static and dynamic spatial inhomogeneities . Nonisovalent chemical substitution in insulating cuprates like $`La_{2x}Sr_xCuO_4`$ and $`YBa_2Cu_3O_{6+x}`$ results in an increase of the energy of the parent phase and creates proper conditions for its competing with other, possibly metallic phases capable to provide an effective screening of the charge inhomogeneity potential.

At the beginning (nucleation regime) a new phase appears in form of a somewhat like a metallic droplet in insulating matrix. The parent phase and droplets of new phase may coexist in the phase separation regime. At different stages of the ”chemical” phase separation regime we deal with the isolated droplets, with the percolation effect and finally with a complete removal of the parent phase . This viewpoint on the phase separation phenomena in cuprates is generally compatible with the pioneer ideas by Emery and Kivelson  and some other model approaches. So, on the basis of the neutron scattering data Egami  conjectured an appearance of a nano-scale heterogeneous structure which is composed of a plenty mobile-carrier existing region of metallic conductivity and semi-localized scarce carrier region with antiferromagnetic spin ordering.

At present, there are numerous experimental indications on the phase separation and nanoscopic stripe-like textures in doped cuprates. Hence when considering the phenomena with a relevant characteristic length, e.g. optical absorption in infrared ($`\lambda 0.52\mu m`$), doped cuprates can be roughly regarded as a binary granular medium, composed of insulating and conducting regions and described in frames of the formalism of an effective medium (EM) theory.

## 2 Physical properties of composites: effective medium approach

In what concerns the optical properties the classical problem of EM is known since long and is related to the calculation of some property of a composite system (typically, the dielectric permittivity $`\epsilon _{eff}`$), those of pure components and their volume fractions been given. Few versions of the theory have been proposed . Many of them were shown to follow from one integral equation  and have inherent shortcoming, assuming the concentration of one of components to be small. Among them the Maxwell-Garnett theory is one of the most popular ones for its simplicity. More sophisticated version of the theory was given by Ping Sheng ; it has been proven valid in entire range of concentration of components.

Metal clusters and small metallic particles embedded in a dielectric matrix have been extensively studied in the past years. One of the most interesting features of these systems is an appearance of the morphological (geometrical), or Mie resonances  related to the excitation of surface plasmon modes in nanoparticle when it scatters the light wave of appropriate frequency. The dispersion law of surface eigenmodes in particles explicitly depends on its geometric shape. Rigorous theory of this effect was first derived by Mie for the case of a spherical particle . In principle, infinite number of geometric resonances should be observed, that correspond to dipolar, quadrupolar etc. surface eigenmodes. However, only dipole contribution is practically accounted. This is valid provided that the particle is placed in homogeneous field, i.e. that the size of the particle is small compared with the wavelength. The frequency of geometric resonance is then easily obtained as the one at which the polarizability of small particle diverges. For the case of spherical particle embedded in the matrix with dielectric permittivity $`\epsilon _d`$ this leads to the equation:

$$\epsilon (\omega )_{part.}+\mathrm{\hspace{0.17em}2}\epsilon _d=\mathrm{\hspace{0.17em}0},$$
(1)
whence the resonance frequency is

$$\omega _r=\frac{\omega _p}{\sqrt{1+\mathrm{\hspace{0.17em}2}\epsilon _d}},$$
(2)
if the Drude’s expression is adopted for the particle permittivity and $`\epsilon _d`$ is assumed to be constant. In the case of arbitrary ellipsoid there are three different principal values of polarizability and the last formula generalizes to

$$\omega _r^i=\omega _p\sqrt{\frac{L_i}{\epsilon _dL_i(\epsilon _d\mathrm{\hspace{0.17em}1})}},i=\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}3},$$
(3)
where $`L_i`$ are three (shape-dependent) depolarizaton factors. In both expressions $`\omega _p`$ is the plasma frequency.

Below we outline in brief some underlying ideas of the EM theory having in mind its further application to doped cuprates. The medium is assumed to be composed of ellipsoidal particles of two kinds. Recalling the lamellar structure of copper oxides we assumed the particles to be randomly oriented in ab plane, one of their principal axes being aligned in c direction. When one addresses the microstructure of such composite two cases may occur, when either of components form the grain and other - the coating. The essential of the theory consists in treating the structure of material as a weighted average of these two conformations in proportion, that depends on volume fraction of components. This has an advantage in correctly displaying the percolation threshold. Final expression which determines an effective dielectric function has a form :

$$fD(\epsilon _{eff},\epsilon _1,\epsilon _2,p)+(1f)D(\epsilon _{eff},\epsilon _2,\epsilon _1,\mathrm{\hspace{0.17em}1}p)=0.$$
(4)
Here, the quantity $`D(\epsilon _{eff},\epsilon _1,\epsilon _2,p)`$ is the orientationally averaged dipole moment of a kind-I particle coated by kind-II particles, embedded in effective medium, $`f`$ and $`(1f)`$ are relative probabilities of occurrence for two kinds of microstructure of the composite. Rough estimate yields $`f`$ in the form: $`f=v_1/(v1+v2)`$, $`v_1=(1\sqrt[3]{p})^3`$, $`v_2=(1\sqrt[3]{1p})^3`$, where $`p`$ is the volume fraction of component I. Analytical derivation of corresponding expression needs the grain and coating to be confocal and leads to the result :

$`D(\epsilon _{eff},\epsilon _1,\epsilon _2,p)`$ $`=`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2}}{}}{\displaystyle \frac{p(\epsilon _1\epsilon _2)\epsilon _2+\left(\left(L_i^{in}pL_i^{out}\right)(\epsilon _1\epsilon _2)+\epsilon _2\right)(\epsilon _2\epsilon _{eff})}{pL_i^{in}(\epsilon _1\epsilon _2)\epsilon _2+\left(\left(L_i^{in}pL_i^{out}\right)(\epsilon _1\epsilon _2)+\epsilon _2\right)\left(\epsilon _{eff}L_i^{out}(\epsilon _{eff}\epsilon _2)\right)}},`$ (5)
where $`L_i^{in,out}`$ are two in-plane depolarization factors of inner and outer ellipsoids, respectively. These are defined as follows:

$$L_i=\frac{a_1a_2a_3}{2}_0^{\mathrm{}}\frac{dt}{(t+a_i^2)\sqrt{_{k=1}^3(t+a_k^2)}},$$
(6)
and can be easily re-expressed in terms of ratios between the semi-axes of ellipsoid.

Despite the shape anisotropy of nanoscopic inclusions the dielectric fuctions for both media I and II together with the effective medium are assumed to be isotropic for a sake of simplicity. This assumption of optically isotropic media implies only a semi-quantitative description of such an anisotropic system as quasi-2D cuprate.

## 3 Optical spectroscopy of doped cuprates: experiment and theoretical models

Optical properties of HTSC cuprates have been the subject of numerous experimental and theoretical investigations. One of the first systematic studies of the doping dependence of the optical reflectance and transmittance for epitaxial films of $`La_{2x}Sr_xCuO_4`$ deposited on $`SrTiO_3`$ substrate were reported in . The most dramatic change with increasing $`x`$ was observed near $`1.5`$ eV. The author speculated that the emergence and near proportional strengthening of the additional absorption centered at $`1.5`$ eV (see Fig. 4) might be related to a dramatic change in the valence band structure. However, the careful measurements of optical reflectivity spectra for bulk single crystals of $`La_{2x}Sr_xCuO_4`$ with $`x`$ varying in wide range, which covers the whole phase diagram of the material, performed by Uchida et al. , reveal only a relatively weak feature near $`1.5`$ eV which becomes apparent only for $`x>0.10`$. The authors suppose that the $`1.5`$ eV feature is extrinsic in origin, and may be related, for instance, to some undetectable amount of oxygen vacancies.

Main effect of doping in $`La_{2x}Sr_xCuO_4`$ and other cuprates is associated with the gradual shift of the spectral weight from the high-energy part of the spectrum to the low-energy part with formation of the prominent midinfrared (MIR) region absorption band and pronounced Drude-like peak. This spectral shift is particularly emphasized by the presence of the isosbestic point at which $`\sigma (\omega )`$ is invariant against $`x`$ .

The authors  state that simple band theory cannot describe the low-energy part of the spectrum within $`CuO_2`$ planes, whereas it nicely mimics the higher-energy excitations. For $`x>0.05`$ the former is clearly separated into a Drude-type peak at $`\omega =0`$ and a broad band centered in the mid-IR region.

One of the puzzling features of the reflectivity spectra in $`La_{2x}Sr_xCuO_4`$ and other cuprates is that the position of the reflectivity edge is almost unchanged against dopant concentration and eventually becomes a plasma edge of the Drude carriers in the overdoped region. In frames of a conventional approach it means that the edge frequency, or a screened plasma frequency, has nearly constant magnitude that rules out the simple relation $`nx`$, which means that the edge is associated straightforwardly with the doped holes. Instead, it is plausible to suppose that $`n(1x)`$ . Comparative analysis of the reflectivity spectra of various cuprates  shows that the plasma frequency of the materials with a nearly half-filled band is almost independent of $`x`$ for small doping but dependent on the average spacing of the $`CuO_2`$ planes. This conclusion clearly contradicts to the results based on the $`tJ`$ model  which show a rather strong dependence of plasma frequency $`\omega _p`$ on hole concentration, especially at small doping.

Extremely lightly oxygen doped crystals $`La_2CuO_{4+y}`$ grown by two different methods (top-seeded and float-zone)  exhibit different low-temperature optical absorption spectra from $`0.1`$ to $`1.5`$ eV photon energy, albeit at room temperature the spectra are nearly identical both each other and those for slightly ”chemically” doped $`La_{2x}Sr_xCuO_4`$  and ”photo-doped” $`La_2CuO_4`$ . In other words, in all cases, irrespective of the ”external” cause, we deal with the same intrinsic mechanism of optical absorption.

Electron-hole optical excitations for the doped cuprate have been calculated by Wagner et al.  via exact diagonalization on finite clusters for a 2D multiorbital Hubbard model. Doping with one and two holes for 12-site lattices introduces a shift of absorption weight from the high-energy side to the low-energy one, in general agreement with experiment. Three additional features can clearly be identified on the low-energy side: first, the Drude peak with weight $`x`$, second, a small feature caused by internal transitions in the Zhang-Rice singlet manifold which may be related to the experimentally observed midinfrared absorption. Finally, a strong absorption structure caused by the transitions from the Zhang-Rice singlet states into the upper Hubbard band was detected and assigned to experimental spectral weight found in the conductivity measurements between the Drude and the charge transfer peaks.

The optical conductivity has been calculated using the Hubbard and $`tJ`$ models by exact diagonalization on small clusters . These studies predict that at low hole density, in addition to a Drude peak with a spectral weight proportional to doping there is a MIR band centered at photon energy $`2J`$, which should be rather broad and presumably caused by hole coupling to spin excitations. The high energy ($`4J`$) of the MIR bands as observed experimentally is in clear contradiction with such a theory.

Later on Lorenzana and Yu  calculated the optical conductivity of $`La_{2x}Sr_xCuO_4`$ in the inhomogeneous Hartree-Fock plus random-phase approximation using the $`pd`$ model with parameters taken from first-principle calculations. The same as for Wagner et al.  CT band peaks now at $`1.7t_{pd}2.7`$ eV. The precursor of the MIR band is dentified with transitions from the localized states on $`Cu`$ to the extended states above the Fermi level. Despite the authors  state close agreement with experiments regarding the peak positions and relative intensities for zero and small doping, their model calculations fail to properly describe the overall evolution of spectra with dielectric-to-metal transition.

## 4 The effective-medium theory of optical conductivity in doped cuprates

We have applied the formalism of effective medium theory  in order to reproduce some tendencies in evolution of optical spectra of high-$`T_c`$ cuprates upon doping. One should note that the optical conductivity and its doping dependence are considered to be one of the most important and informative characteristics of the doped cuprates.

In particular, we have attempted to describe the experimental data on doping dependence of optical conductivity in $`La_{2x}Sr_xCuO_4`$  in framework of the EM model. It should be emphasized that in our case the EM theory is nothing but a phenomenological approach to the description of the complex spatially non-uniform rearrangements in the electronic structure of doped cuprate. The concept of binary composite merely reflects the fact that upon doping some fraction of electrons remains localized in the regions of parent insulating phase while other electrons becomes nearly metallic inside regions embedded into parent insulating matrix.

It is reasonable to conjecture, that different electronic states of La2CuO4 do not contribute equally to this process. In particular, we simply assumed that CT excitation band(s) peaked near $`2.5`$ eV is involved in such a ”metallization” while the states at higher energy remain rigid. These ideas can easily account for the behaviour of optical conductivity of La2-xSrxCuO4: upon doping the spectral weight rapidly transfers from the CT band(s) to a narrow band of extended states, while the high-energy tails converge beyond $`\mathrm{\hspace{0.17em}4.5}`$ eV for different $`x`$.

The imaginary part of $`\epsilon `$ for the parent undoped cuprate $`La_2CuO_4`$ was taken from  and . The spectrum of Im $`\epsilon `$ of $`La_2CuO_4`$ we used in our calculations was modeled by a superposition of two CT bands fitted by a sum of two Gaussians centered at $`2.0`$ and $`3.0`$ eV, and the ”rigid” high-energy contribution fitted by a single Lorentzian centered at $`6.0`$ eV. Only the first CT term is involved in EMT computations, its spectral weight being partly substituted by that of new metallic-like phase. As for the real part of dielectric function, near dispersionless value close to $`ϵ=\mathrm{\hspace{0.17em}4}`$ has been chosen.

To model the dielectric function of the metallic phase we used the simplest Drude form

$$\epsilon (\omega )=\epsilon _{\mathrm{}}\left(1\frac{\omega _p^2}{\omega (\omega +i\gamma )}\right),$$
(7)
where $`\omega _p`$ and $`\gamma `$ are the plasma frequency and the electron scattering rate, respectively. Together with $`\epsilon _{\mathrm{}}`$ these are regarded as adjustable parameters which magnitude does not depend on doping, frequency, and temperature. The results presented below were obtained with the values: $`\omega _p=\mathrm{\hspace{0.17em}1.65}`$ eV, $`\gamma =\mathrm{\hspace{0.17em}0.6}`$ eV, $`\epsilon _{\mathrm{}}=\mathrm{\hspace{0.17em}1.0}`$, which once tuned, were retained unchanged against the doping level.

Another set of free parameters stand for the microtexture of the sample under investigation. Here we take into account numerous manifestations of stripe-like textures in doped cuprates. We denote the in-plane semi-axes ratio for particles of both kinds as $`\alpha `$ and that of out-of-plane and major in-plane semi-axes as $`\beta `$. These geometrical parameters enter the theory through the depolarization factors (6) and were found to affect strongly the dielectric function of doped cuprate.

Calculated spectra of optical conductivity $`\sigma (\omega )`$ for different level of doping $`x`$ are shown in Fig. 1 together with experimental curves from . The numbers near the curves stand for the fractions of metal phase $`p`$, which are functions of $`x`$. The geometrical parameters $`\alpha _1=\mathrm{\hspace{0.17em}0.3},\beta _1=\mathrm{\hspace{0.17em}0.09}`$ and $`\alpha _2=\mathrm{\hspace{0.17em}1.0},\beta _2=\mathrm{\hspace{0.17em}0.05}`$ were adopted for metallic and insulating grains, respectively. This implies the effective metallic nanoparticle to be stripe-shaped (elongated in ab \- plane and nearly plate in the c \- direction) and dielectric ones to be oblate ellipsoid-like round pellets (equal in-plane semi-axes and much smaller out-of-plane semi-axis). Spectral weight was found to transfer dramatically from CT band to a low-energy band at $`0.7`$ eV in the range $`p\mathrm{\hspace{0.17em}0.3}`$ that corresponds to $`x\mathrm{\hspace{0.17em}0.06}`$. Further increase of metal fraction results in full removal of $`2.2`$ eV feature of parent compound and emergence of a new weak peak at $`\mathrm{\hspace{0.17em}1.5}`$ eV for large enough $`p`$, really seen in the heavily doped compositions ($`x\mathrm{\hspace{0.17em}0.15}`$). Simultaneously, the $`0.7`$ eV band shifts to lower energies and gradually transforms into the Drude-like peak. Many doped cuprates demonstrate isosbestic behaviour of optical conductivity, particularly well pronounced e.g. in $`T^{}`$ phase of $`Nd_{2x}Ce_xCuO_4`$: $`\sigma (\omega ,x)`$ cross the same point for various $`x`$. It is seen that our model theory remarkably well reproduces this feature.

Overall agreement of simple EM theory with experimental data  is rather impressive despite such oversimplified assumptions as invariable value of shape and size of grains, plasma frequency, relaxation rate and other quantities regardless the doping level. At the same time one observe clear shortcomings in what concerns the detailed quantitative description of the low-frequency ($`\mathrm{}\omega \gamma `$) spectral range and the percolation phenomenon. The EM theory yields the percolation threshold volume fraction $`p_c0.6`$ that seems too high. The more detailed description of the low-frequency spectral range in frames of the EM theory first of all implies the account of the size and shape distribution of metallic/dielectric grains along with the distribution of quantities like $`\omega _p`$ and $`\gamma `$.

We have identified the low-energy peak and $`1.5`$ eV feature with geometric Mie resonances as they coincide with divergencies in some terms of polarizability of composite (nullification of denominator in (5)). This is a quite novel type of resonant absorption, inherent for composite materials, which is not originated straightforwardly from either electronic transitions. The strong influence of the grain shape on the $`1.5`$ eV feature is fairly well illustrated in Fig. 2 which demonstrates the calculated optical conductivity $`\sigma (\omega )`$ given the metallic volume fraction $`p=0.55`$ at different values of parameter $`\beta `$ which defines the ratio of the out-of-plane and in-plane semi-axes of dielectric grain.

The prominent $`\mathrm{\hspace{0.17em}1.5}`$ eV peak is substantially more pronounced in optical absorption spectra reported by Suzuki  for thin-film samples. We show that good agreement with both sets of the experimentals  could be obtained merely adjusting the ”geometric”, or shape parameters of metallic and insulating droplets without modifying any underlying physics of the model. Calculated spectra of optical conductivity $`\sigma (\omega )`$, absorption $`\alpha (\lambda )`$ and transmittance $`T(\lambda )`$ for the $`0.6\mu `$m thin film samples of $`La_{2x}Sr_xCuO_4`$ are presented in Fig. 3 and Fig. 4. Actually, all physical parameters were retained unvaried except of Drude’s relaxation rate: $`\gamma =\mathrm{\hspace{0.17em}0.4}`$ eV was set instead of previous value 0.6 eV, but new values for shape parameters $`\alpha _1=\mathrm{\hspace{0.17em}0.9},\beta _1=\mathrm{\hspace{0.17em}0.01}`$ and $`\alpha _2=\mathrm{\hspace{0.17em}0.9},\beta _2=\mathrm{\hspace{0.17em}1.5}`$ are taken. This implies the oblate ellipsoid-like metallic grains and predominantly 3D character of dielectric regions. The most essential consequence of these changes is associated with pronounced development of $`1.5`$ eV feature in optical conductivity. Accordingly, the $`\alpha (\lambda )`$ and $`T(\lambda )`$ spectra fit fairly well the experimental curves .

Real part of optical conductivity addressed above relates to imaginary part of dielectric permittivity. As an important additional check of validity of our simple model we address the calculation of the electron energy dependent loss function Im$`(1/ϵ(\omega ))`$ which is determined both by real and imaginary parts of $`ϵ(\omega )`$. The Figure 5 presents a calculated spectral dependence of the electron energy loss function for a model $`La_{2x}Sr_xCuO_4`$ system given the same model parameters used above for description of the optical conductivity for the bulk samples (see Fig. 1). Comparison with experimental data  shows that our simple model allows correctly describe all features of the EELS spectra, in particular, a non-monotonous doping dependence of the peak in the loss function. It should be noted that in contrast with conventional approach both reflectivity edge and the main peak in the loss function in frames of EM theory do not relate straightforwardly to the bare plasma frequency $`\omega _p`$.

We assumed through the article that the doping level $`x`$ (say, strontium content) can be simulated by an appropriate tuning of metal volume fraction $`p`$ in the composite. It is, however, desirable to establish the relation between these quantities in more accurate manner. For the case of optical conductivity in bulk single-crystalline $`La_{2x}Sr_xCuO_4`$ samples  this can be done e.g. through the calculation relative decrease of the $`2.2`$ eV peak with doping, then adjusting $`p`$ to keep the same proportion in theoretical graph. The Fig. 6 illustrates the relation $`p(x)`$ obtained in this way. The saturation behaviour, discernible for large enough $`x`$ can be attributed to the percolation phenomenon.

## 5 Conclusion

In conclusion, we have proposed a new approach to theoretical description of doped cuprates, treating these as inhomogeneous composite materials, containing the dielectric and metallic stripe-like nanoparticles. The formalism of effective medium theory is then applied for calculation of dielectric permittivity and optical spectra of $`La_{2x}Sr_xCuO_4`$ with $`x`$ varying in a wide range. Reasonable quantitative agreement with experiment has been obtained. The model in its simplest form was found able to reproduce all essential features of the transmittance , optical conductivity $`\sigma (\omega )`$, and EELS spectra . Substantial difference in spectral and doping dependence of optical conductivity for thin-film  and bulk samples  is easily explained if to only assume different shape of metallic and dielectric regions in both materials. New peaks in $`\sigma (\omega )`$ and absorption spectra, that emerge in the midinfrared range upon doping are attributed to geometrical (Mie’s) resonances. Overall, the model theory allows to properly describe the evolution of optical spectra at the percolative insulator-to-metal transition.

Finally, we would like to note that an occurrence of the phase separation with percolative nature of the insulator-to-metal transition is accompanied by many specific features which in many cases mask the real electronic structure and can lead to the erroneous conclusions. In particular, making use of either experimental data as a trump in favor of either mechanism of the high-$`T_c`$ superconductivity should be made with some caution, if a phase homogeneity of the samples under examination is questionable. In this connection one may mention a long-standing discussion regarding the nature of the midinfrared bands in doped cuprates.

Further extension of the model developed here should include the size and shape distribution of metallic-like inclusions, temperature effects, as well as effects of different external factors such as pressure, isotopic substitution, photo-doping. Of course, this semi-empirical model needs in more detailed microscopic reasoning.

We are grateful to N. Loshkareva, Yu. Sukhorukov, K. Kugel, P. Horsch for discussions. One of the authors (A.S.M.) acknowledges stimulating discussion with N. DelFatti. The research described in this publication was made possible in part by Award No.REC-005 of the U.S. Civilian Research & Development Foundation for the Independent States of the Former Soviet Union (CRDF). The authors acknowledge a partial support from the Russian Ministry of Education, grant E00-3.4-280.


text:
# Modulational instability in nonlocal nonlinear Kerr media

## I Introduction

Modulational instability (MI) constitutes one of the most fundamental effects associated with wave propagation in nonlinear media. It signifies the exponential growth of a weak perturbation of the amplitude of the wave as it propagates. The gain leads to amplification of sidebands, which breaks up the otherwise uniform wave front and generates fine localized structures (filamentation). Thus it may act as a precursor for the formation of bright spatial solitons. Conversely the generation of dark spatial solitons requires the absence of MI of the constant intensity background.

The phenomena of MI has been identified and studied in various physical systems, such as fluids , plasma , nonlinear optics , discrete nonlinear systems, such as molecular chains  and Fermi-resonant interfaces and waveguide arrays , etc. It has been shown that MI is strongly affected by various mechanisms present in nonlinear systems, such as higher order dispersive terms in the case of optical pulses , saturation of the nonlinearity , and coherence properties of optical beams .

In this work we study the MI of plane waves propagating in a nonlinear Kerr type medium with a nonlinearity (the refractive index change $`\mathrm{\Delta }n`$, in nonlinear optics) that is a nonlocal function of the incident field. We consider a phenomenological model, in which the nonlocal nonlinearity, induced by a wave (e.g. an optical beam) with the intensity $`I(x,z)`$, can be represented in general form as

$$\mathrm{\Delta }n(I)=s_{\mathrm{}}^{\mathrm{}}R(x^{}x)I(x^{},z)𝑑x^{},$$
(1)
where $`x`$ is the transverse spatial coordinate and $`s`$=1 ($`s`$=$`1`$) corresponds to a focusing (defocusing) nonlinearity. The evolution coordinate $`z`$ can be time, as for Bose-Einstein Condensates (BEC’s), or the spatial propagation coordinate, as for optical beams. We consider only symmetric spatial response functions that are positive definite and (without loss of generality) obey the normalization condition

$$_{\mathrm{}}^{\mathrm{}}R(x)𝑑x=1.$$
(2)
Thus we exclude asymmetric effects, such as those generated by asymmetric temporal response functions (with $`x`$ being time), as in the case of the Raman effect on optical pulses .

In nonlinear optics Eq. (1) represents a general phenomenological model for media in which the nonlinear refractive index change (or polarization) induced by an optical beam is determined by some kind of a transport process. It may include, e.g., heat conduction in materials with a thermal nonlinearity  or diffusion of molecules or atoms accompanying nonlinear light propagation in atomic vapours . Nonlocality also accompanies the propagation of waves in plasma , and a nonlocal response in the form (1) appears naturally as a result of many body interaction processes in the description of Bose-Einstein condensates .

The width of the response function $`R(x)`$ relative to the width of the intensity profile $`I(x,z)`$ determines the degree of nonlocality, as illustrated in Fig. 1. In the limit of a singular response, $`R(x)`$=$`\delta (x)`$ \[see Fig. 1(a)\], the nonlinearity is a local function of the intensity,

$$\mathrm{\Delta }n=sI,$$
(3)
i.e. the refractive index change at a given point is solely determined by the wave amplitude at that very point. For this well-known Kerr nonlinearity, which appears in all areas of physics, MI depends only on the sign $`s`$ .

With increasing width of $`R(x)`$ the intensity in the vicinity of the point $`x`$ also contributes to the index change at that point. For a weak nonlocality, when the width of $`R(x)`$ is finite but still small compared to the width of $`I(x,z)`$ \[see Fig. 1(b)\], one can expand $`I(x^{},z)`$ around $`x^{}`$=$`x`$ in (1) and obtain the simplified model

$$\mathrm{\Delta }n=s(I+\gamma _x^2I),\gamma =\frac{1}{2}_{\mathrm{}}^{\mathrm{}}x^2R(x)𝑑x,$$
(4)
where $`\gamma `$ is a small positive parameter. This diffusion type model of the nonlocal nonlinearity is a model in its own right in plasma physics, where $`\gamma `$ can take any sign . It was also applied to BEC’s , nonlinear optics  and in the continuum limit of the theory of energy transfer in biomolecular systems . In contrast to the local Kerr limit the MI now depends not only on the sign $`s`$, but also on the intensity of the plane wave .

In the limit of a strongly nonlocal response, where $`R(x)`$ is much broader than the intensity profile \[see Fig. 1(d)\], one can expand the response function around $`x^{}`$=$`x`$ in (1) and obtain the simple linear model

$$\mathrm{\Delta }n=s(c_0+c_1x+c_2x^2),$$
(5)
where $`c_{02}`$ are constants. Since this model is linear all plane waves are stable, - there is no MI. The evolution of optical beams in such a medium was considered in the case when $`c_1`$=0 , which requires the center of mass of the beam to be always zero.

Modulational instability has thus been studied in the different limits. However, in the general case (1) the problem of MI has only been studied for a few particular cases of the response function $`R(x)`$. Here we present an analytical study of the general case with arbitrary profile $`R(x)`$ and confirm the results for three specific examples. We further show that the MI results of the weakly nonlocal model (4) lead to erroneous conclusions when used to predict the behavior of the full model (1).

## II General Theory

We will consider an optical beam propagating along the $`z`$-axis, with the scalar amplitude of the electric field expressed as

$$E(x,z)=\psi (x,z)\mathrm{exp}(iKzi\mathrm{\Omega }t)+c.c.,$$
(6)
where $`K`$ is the wavenumber, $`\mathrm{\Omega }`$ is the optical frequency and $`\psi (x,z)`$ is the slowly varying amplitude. Substituting Eqs. (1) and (6) into Maxwell’s equations we obtain the nonlocal nonlinear Schrödinger (NLS) equation

$$i_z\psi +\frac{1}{2}_x^2\psi +s\psi _{\mathrm{}}^{\mathrm{}}R(x^{}x)I(x^{},z)𝑑x^{}=0,$$
(7)
where $`I(x,z)`$=$`|\psi (x,z)|^2`$ is the intensity of the beam. The model (7) permits plane wave solutions of the form

$$\psi (x,z)=\sqrt{\rho _0}\mathrm{exp}(ik_0xi\omega _0z),\rho _0>0,$$
(8)
where $`\rho _0`$, $`k_0`$, and $`\omega _0`$ are linked through the nonlinear dispersion relation

$$\omega _0=\frac{1}{2}k_0^2s\rho _0,$$
(9)
which is the same as for the standard local NLS equation

$$i_z\psi +\frac{1}{2}_x^2\psi +s|\psi |^2\psi =0.$$
(10)
Next, let us carry out a linear stability analysis of the plane wave solutions (8). Assume that

$$\psi (x,z)=[\sqrt{\rho _0}+a_1(x,z)]\mathrm{exp}(ik_0xi\omega _0z),$$
(11)
where $`a_1(x,z)`$ is a small complex perturbation. Inserting this expression into the nonlocal NLS equation (7) and linearizing around the solution (8) yields the evolution equation for the perturbation

$`i_\tau a_1+{\displaystyle \frac{1}{2}}_\xi ^2a_1`$ (12)
$`+2s\rho _0{\displaystyle _{\mathrm{}}^{\mathrm{}}}R(\xi ^{}\xi )\mathrm{Re}\{a_1(\xi ^{},\tau )\}𝑑x^{}=0.`$ (13)
In deriving this equation we have used the nonlinear dispersion relation (9) and gone into a coordinate frame moving with the group velocity $`c_g=d\omega _0/dk_0=k_0`$, i.e.

$$\tau =z,\xi =xc_gz$$
(14)
Decomposing the perturbation into real and imaginary parts, $`a_1=u+iv`$, we obtain two coupled equations

$$\begin{array}{c}_\tau u+\frac{1}{2}_\xi ^2v=0,\\ \\ _\tau v\frac{1}{2}_\xi ^2u2s\rho _0_{\mathrm{}}^{\mathrm{}}R(\xi ^{}\xi )u(\xi ^{},\tau )𝑑\xi ^{}=0.\end{array}$$
(15)
By introducing the Fourier transforms

$`\widehat{u}(k,\tau )`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}u(\xi ,\tau )\mathrm{exp}[ik\xi ]𝑑\xi ,`$ (16)
$`\widehat{v}(k,\tau )`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}v(\xi ,\tau )\mathrm{exp}[ik\xi ]𝑑\xi ,`$ (17)
$`\widehat{R}(k)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}R(\xi )\mathrm{exp}[ik\xi ]𝑑\xi ,`$ (18)
and exploiting the convolution theorem for Fourier transforms, the linearized system is converted to a set of ordinary differential equations in $`k`$-space

$$\begin{array}{c}_\tau \widehat{u}\frac{1}{2}k^2\widehat{v}=0,\\ \\ _\tau \widehat{v}+\frac{1}{2}k^2\widehat{u}2s\rho _0\widehat{R}\widehat{u}=0,\end{array}$$
(19)
which can be written in the compact matrix form

$$_\tau \underset{¯}{X}=\underset{¯}{\underset{¯}{A}}\underset{¯}{X},$$
(20)
where the vector $`\underset{¯}{X}`$ and matrix $`\underset{¯}{\underset{¯}{A}}`$ are defined as

$$\underset{¯}{X}=\left[\begin{array}{c}\widehat{u}\\ \widehat{v}\end{array}\right],\underset{¯}{\underset{¯}{A}}=\left[\begin{array}{cc}0& \frac{1}{2}k^2\\ \frac{1}{2}k^2+2s\rho _0\widehat{R}(k)& 0\end{array}\right].$$
(21)
The eigenvalues $`\lambda `$ of the matrix $`\underset{¯}{\underset{¯}{A}}`$ are given by

$$\lambda ^2=k^2\rho _0\left[\alpha k^2s\widehat{R}(k)\right],$$
(22)
where we have defined the parameter $`\alpha `$ as

$$\alpha =1/(4\rho _0)$$
(23)
Notice that if the response function $`R(x)`$ is real and symmetric then so is the Fourier spectrum of $`R(x)`$, i.e. $`\widehat{R}(k)=\widehat{R}^{}(k)=\widehat{R}(k)`$.

The general eigenvalue equation (22) constitutes the basis of our further study. First of all we notice that $`\widehat{R}(0)`$=1, since the response function is normalized according to Eq. (2). For a focusing nonlinearity ($`s`$=$`+1`$) we therefore always have, by continuity, that $`\alpha k^2s\widehat{R}(k)<0`$, and thus $`\lambda ^2>0`$, in a certain wavenumber band symmetrically centered about the origin, where $`k`$ is sufficiently small. It follows that we always will have (long wave) MI in the focusing case, independently of the details in the behavior of the response function.

In contrast, in the defocusing case ($`s`$=$`1`$), the stability properties depend in a sensitive way on the response function. Thus the nonlocality could possibly lead to MI, which would never be present in the defocusing local NLS equation (10).

The well-known modulational instability (stability) result for the standard local NLS equation (10) is easily recovered from the general eigenvalue equation (22) by setting $`R(x)`$=$`\delta (x)`$, where $`\delta (x)`$ is the Dirac delta function: We get

$$\lambda ^2=k^2\rho _0(\alpha k^2s)$$
(24)
where $`s=+1`$ $`(s=1)`$ yields instability (stability).

## III Weakly nonlocal limit

The MI properties in the local limit (3), described by the NLS equation (10), are well-known and the strongly nonlocal limit (5) is linear, so obviously it does not display MI. Here we briefly consider the interesting regime of a weak nonlocality, i.e., the situation when the typical width of the response kernel is small compared to the characteristic length of the wave modulation (see Fig. 1). In this case the nonlinearity has the form (4) and the NLS equation (7) reads

$$i_z\psi +\frac{1}{2}_x^2\psi +s(|\psi |^2+\gamma _x^2|\psi |^2)\psi =0,$$
(25)
where the nonlocality or diffusion parameter $`\gamma `$ is defined in Eq. (4). Equation (25) has been discussed in the literature in the context of plasma physics  and the continuum limit of discrete molecular structures  and has been shown to posses bright and dark soliton solutions. Their exact analytical soliton form was found recently . Importantly, $`\gamma `$ is small and positive when Eq. (25) represents the weakly nonlocal limit of the general nonlocal model (7), whereas $`\gamma `$ can take any sign in plasma physics, where Eq. (25) is a model in its own right.

The instability condition for the weak nonlocality is most easily obtained by Taylor expanding the spectrum of the response kernel about $`k`$=0 and only retain the two lowest order nontrivial terms in the expansion, i.e.

$$\widehat{R}(k)1\gamma k^2,$$
(26)
where obviously $`\gamma k^21`$ must be fulfilled. The eigenvalue equation then reads

$$\lambda ^2=k^2\rho _0[(\alpha +s\gamma )k^2s].$$
(27)
This result does not depend on the detailed behavior of the response kernel, as should be expected.

The eigenvalue Eq. (27) clearly shows that MI is always present for a focusing nonlinearity ($`s`$=$`+1`$), since for sufficiently small wavenumbers $`\lambda ^2>0`$ for all values of $`\alpha `$ and $`\gamma `$. This is in correspondance with the results of the general nonlocal model (7).

More interesting is the defocusing case ($`s`$=$`1`$), where there is a critical value $`\rho _0^{\mathrm{cr}}`$ of the intensity,

$$\rho _0^{\mathrm{cr}}=1/(4\gamma ),$$
(28)
which depends on the degree of nonlocality through the parameter $`\gamma `$, and corresponds to $`\alpha =\gamma `$. When $`\gamma `$ is positive (and $`s=1`$) low intensity plane waves with $`\rho _0\rho _0^{\mathrm{cr}}`$ are modulationally stable, whereas MI appears at high intensities, $`\rho _0>\rho _0^{\mathrm{cr}}`$, despite the nonlinearity being defocusing. This MI is short-wave, appearing for pertubations with wavenumbers exceeding the threshold value

$$k^{\mathrm{cr}}=1/\sqrt{|\alpha \gamma |}.$$
(29)
The stability criterion $`\rho _0<\rho _0^{\mathrm{cr}}`$ corresponds exactly to the criterion for existence of dark soliton solutions in this case . The existence of MI in the weakly nonlocal model is in sharp contrast to the local defocusing Kerr nonlinearity, for which there is no MI.

In the case $`\gamma <0`$ in Eq.(25), relevant in plasma physics, nonlocality does not affect the MI properties: For a focusing (defocusing) nonlinearity all plane waves are modulationally unstable (stable). In Table. I we have summarized the MI results of the weak nonlocality model.

So only for defocusing media with $`\gamma >0`$ does a weak nonlocality change the MI properties in terms of existence, as compared to the local NLS limit. However, nonlocality does always affect the structure of the MI gain spectrum, as illustrated in Fig. 2 for $`\gamma >0`$. In the focusing case Fig. 2 clearly shows how the nonlocality tends to suppress MI by decreasing both the maximum gain and the gain bandwidth. In contrast, in the defocusing case, the nonlocality promotes MI, - the higher the degree of nonlocality $`\gamma `$ the lower the threshold value of the wavenumber, above which the instability develops.

In Fig. 3 we show the development of MI of a plane wave in a self-defocusing weakly nonlocal medium, as obtained by numerical integration of the NLS equation (25). As an initial condition we used a plane wave with periodically perturbed amplitude. It is evident that this perturbation grows as the wave propagates.

Although Eqs.(28) and (29) predict MI in the defocusing case it is important to stress that this result does not actually apply to the case of weak nonlocality. In deriving the model (25) and the expansion (26) it was assumed that

$$\gamma k^21,$$
(30)
which means that the width of the response function must be small compared to the characteristic length-scale of the modulation. However, from Eq. (27) we get that MI occurs when

$$(\alpha \gamma )k^2+1<0,$$
(31)
which can be rewritten as

$$\gamma k^2>1+\alpha k^2.$$
(32)
Since $`\alpha =1/(4\rho _0)`$ is positive this condition cannot be fulfilled without violating the perturbation requirement (30). Hence, in the region of parameters where MI occurs (and dark solitons do not exist), the weak nonlocality approximation breaks down. This is in accordance with the analysis of the full nonlocal model (7), which, as we show below, does not display MI for self-defocusing nonlinearities with a smooth nonlocal response function.

It should be stressed however that the appearance of MI is a valid result in the situation when the weakly nonlocal model (7) is a model in its own right, as for example for upper hybrid modes in plasma .

## IV Fully nonlocal response - examples

We will now discuss the stability issue in the general nonlocal Kerr nonlinearity with an arbitrary degree of nonlocality by considering a few specific examples of the nonlocal response function.

### A Gaussian response function

Let the response function of the nonlocal nonlinear medium have the following Gaussian form

$$R(x)=\frac{1}{\sigma \sqrt{\pi }}\mathrm{exp}\left[\frac{x^2}{\sigma ^2}\right].$$
(33)
In this case $`\gamma =\sigma ^2/2`$, the Fourier transform of the response function is given by

$$\widehat{R}(k)=\mathrm{exp}\left[\frac{1}{4}\sigma ^2k^2\right],$$
(34)
and the eigenvalue equation (22) reduces to

$$\lambda ^2=k^2\rho _0\left(\alpha k^2s\mathrm{exp}\left[\frac{1}{4}\sigma ^2k^2\right]\right).$$
(35)
For the defocusing case ($`s=1`$) we find modulational stability for all values of the intensity $`\rho _0`$, in contrast to the behaviour predicted by the weakly nonlocal model (25), according to which plane waves become unstable if their intensity is higher than the critical intensity $`\rho _0^{\mathrm{cr}}`$. Thus, as expected, this proves that the weakly nonlocal model gives incorrect predictions of the MI properties of the full general nonlocal model.

A simple analysis reveals finite band width MI in the focusing case ($`s=+1`$) with a growth rate given by

$$\left|\text{Re}\lambda \right|=\left|k\right|\sqrt{\rho _0}\sqrt{\mathrm{exp}\left[\frac{1}{4}\sigma ^2k^2\right]\alpha k^2}$$
(36)
The dependence of the growth rate on the wave number and degree of nonlocality $`\sigma `$ is displayed in Fig. 4. Evidently, as the degree of the nonlocality increases the MI bandwidth shrinks and the maximum growth rate decreases. Thus nonlocality clearly tends to suppress MI in this case, as also predicted by the weakly nonlocal model.

In Fig. 5 we demonstrate the influence of nonlocality on the dynamic development of MI. To this end we integrated numerically the NLS equation (7) with the Gaussian response function (33). As an initial conditions we used a unit intensity plane wave with imposed weak periodic perturbation

$$\psi (x,z=0)=1+10^4\mathrm{cos}(1.5x)$$
(37)
It is evident that increasing the degree of nonlocality $`\sigma `$ suppresses the instability by lowering the gain.

### B Exponential response function

Another interesting and important example involves the exponential response function

$$R(x)=\frac{1}{2\sigma }\mathrm{exp}\left[\frac{|x|}{\sigma }\right],$$
(38)
for which $`\gamma =2\sigma ^2`$, with $`\sigma `$ again determining the degree of nonlocality. The importance of this particular response lies in the fact that it applies exactly to media with a nonlinearity (i.e. nonlinear refractive index $`\mathrm{\Delta }n(x)`$) described by the diffusion equation

$$\sigma ^2\frac{^2\mathrm{\Delta }n}{x^2}+\mathrm{\Delta }n=|\psi |^2.$$
(39)
The equivalence between the nonlocal nonlinearity (1) with the exponential response (38) and the diffusion model (39) is easily shown by Fourier transformation.

As discussed by Litvak et al.  the model (39) describes the nonlinear thermal self-focusing of electromagnetic waves in a weakly ionized plasma. In this particular case $`\mathrm{\Delta }n`$ corresponds to the change of the plasma temperature. From the Fourier transform

$$\widehat{R}(k)=\frac{1}{1+\sigma ^2k^2},$$
(40)
we obtain the eigenvalue equation

$$\lambda ^2=k^2\rho _0\left[\alpha k^2s\frac{1}{1+(\sigma k)^2}\right],$$
(41)
which coincides with that derived in Ref.  \[Eq. (4.1) in this reference\]. The dispersion relation (41) imposes MI properties similar to those found for the Gaussian response. In particular, it follows that in the self-defocusing regime ($`s=1`$) the plane wave solution is always modulationally stable. On the other hand, in the self-focusing case plane wave develops transverse instability with the growth rate given by

$$|\text{Re}\lambda |=|k|\sqrt{\rho _0}\sqrt{\frac{1}{1+(\sigma k)^2}\alpha k^2}.$$
(42)
In Fig. (6) we show the MI growth rate versus the wave number for a few values of the nonlocality parameter $`\sigma `$. As for the Gaussian response the nonlocality tends to suppress MI by narrowing the gain bandwidth and decreasing the maximum gain. However, only in the extreme limit of an infinitely nonlocal response ($`\sigma \mathrm{}`$) does MI completely disappear, as the medium no longer exhibits any nonlinearity.

### C Rectangular response function

Finally we consider the rectangular response function

$$R(x)=\{\begin{array}{ccc}\frac{1}{2\sigma }& \text{for}& |x|\sigma \\ & & \\ 0& \text{for}& |x|>\sigma \end{array}$$
(43)
for which $`\gamma =\sigma ^2/3`$. From the Fourier transform

$$\widehat{R}(k)=\frac{\mathrm{sin}(k\sigma )}{k\sigma },$$
(44)
we obtain the eigenvalue equation

$$\lambda ^2=k^2\rho _0\left[\alpha k^2s\frac{\mathrm{sin}(k\sigma )}{k\sigma }\right].$$
(45)
In the focusing case one again finds that there is always (for any value of $`\alpha `$ and $`\sigma `$) a finite bandwidth MI gain band symmetrically located about the origin in $`k`$-space, i.e. for $`|k|k_{\mathrm{cr}}`$, where

$$\alpha k_{\mathrm{cr}}^2=\frac{\mathrm{sin}(k_{\mathrm{cr}}\sigma )}{k_{\mathrm{cr}}\sigma },$$
(46)
This is illustrated in Fig. 7, where we show the gain profiles for several degrees of nonlocality $`\sigma `$. Again MI is being suppressed by the nonlocality, but now additional MI bands appear for large values of $`\sigma `$ (see $`\sigma `$=20).

In the defocusing case a careful analysis reveals that for large and moderate values of the parameter $`\alpha /\sigma ^2`$, i.e. small and moderate values of $`\rho _0\sigma ^2`$, plane waves are modulationally stable, while in the opposite regime one or more MI bands appear, as shown in Fig. 8. This instability is caused by the negative-sign bands of the sinc$`(x)`$-profile Fourier Transform (44), with the first instability band being in the region $`\pi k\sigma 2\pi `$. More specifically the instability appears when the background intensity exceeds the threshold value

$$\rho _0^{\mathrm{th}}>21.05/\sigma ^2.$$
(47)
## V Conclusions

We studied the modulational stability properties of plane waves in nonlinear media with a general nonlocal Kerr nonlinearity, described by a nonlocal NLS equation. We derived an analytical expression for the growth rate of the instability for arbitrary response profiles, i.e. for arbitrary degrees of nonlocality. We found that in the self-focusing case the plane wave solution is always unstable, - the nonlocality tends to suppress the instability, independently of the form of the response function, but can never suppress is completely.

In self-defocusing media the stability properties depend strongly on the specific form of the nonlocal response. For any strictly positive continuous response profile all plane waves are stable for any degree of nonlocality.However, exotic profiles with jump-discontinuities like the rectangular response, may actually promote MI when the intensity is sufficiently high”

In the so-called weakly nonlocal limit (narrow response function compared to the intensity profile) we find that plane waves become modulationally unstable if their intensity exceeds a certain critical value. However, we show that this surprising prediction of MI in the defocusing case can never apply to the full nonlocal model with a smooth response profile. Thus one should be cautious with applying results of the simple weakly nonlocal model to the real general case, in particular when studying MI and dark solitons.

This work was supported by the Danish Technical Research Council (STVF - Talent Grant 5600-00-0355), the Danish Natural Sciences Foundation (SNF - grant 9903273), and the Graduate School in Nonlinear Science (The Danish Research Academy).