text:
# The exosystems around HD 169830 and HD 12661: are they dynamical twins?

## 1 Introduction

The radial velocity survey for extrasolar planetary systems conducted by the Geneva Extrasolar Planet Search Team has recently revealed a new planetary system around HD 169830. This discovery has been announced during the XIX-th IAP Colloquium ”Extrasolar Planets: Today & Tomorrow”, (Paris, June 30 - July 4, 2003). In this paper we perform a preliminary analysis of the orbital stability of the HD 169830 system and its dependence on the orbital parameters that are unconstrained by radial velocity (RV) observations. In the numerical experiments, we use two sets of the orbital parameters from the 2-Keplerian fits to the RV measurements. The first set was published by the authors on their web site in June 2003 111http://obswww.unige.ch/~udry/planet/hd169830\_syst.html and the second, most recent set is from Mayor et al. (2003). Our early analysis of the first set of the orbital elements revealed a curious, close similarity to the orbital elements of the HD 12661 system Fischer et al. (2003); Goździewski & Maciejewski (2003). These two hierarchical planetary systems (HPS) appeared to have very similar ratios of the minimal masses $`\mu =m_\text{b}/m_\text{c}1.3`$, the semi-major axes $`\alpha =a_\text{b}/a_\text{c}0.3`$ and similar initial eccentricity of the outer planet $`e_\text{c}0.3`$ with simultaneously small eccentricity of the inner planet $`e_\text{b}0`$. As we show in this paper, in such a case both systems would be located in a large island of the secular apsidal resonance with the semi-major axes antialigned in the exact resonance. There could be another planetary system, around HD 160691 Jones et al. (2002), with a qualitatively similar dynamics as suggested by Goździewski et al. (2003). In a relatively small number of about 14 multi-planetary systems known to the date222http://www.encyclopaedia.fr, these three cases of supposedly identical configuration inspired us to ask whether such a dynamical resemblance can be casual or has a deeper origin in their formation scenario and/or orbital evolution.

However recently, the Geneva team has published a new orbital solution substantially different than the first one (see Mayor et al., 2003). In the new solution the elements of the outer companion are dramatically altered. The period and the semi-major axis of the outer planet are changed by $`600`$ d and $`0.7`$AU respectively. This obviously leads to a qualitatively different orbital configuration than the one analyzed before. The proximity of the HD 169830 and HD 12661 systems does not seem to be so clear anymore. Still, the orbital elements of the inner companions appear to be almost identical in the two planetary systems. Because the observational window of HD 169830 covers only about 1 orbital period of the outer planet, its orbital elements are still not well constrained. In our numerical experiments, we use both solutions although the old one essentially for a reference and an interesting comparison with the predictions given by the secular octupole-level theory (Lee & Peale, 2003).

According to Lee & Peale (2003), one has to employ a proper representation of the orbital fits to the RV observations, particularly in the case of the HPS. The orbital fits to the RV of HPS are best interpreted in Jacobi coordinates and these coordinates should be used in the dynamical investigations. Specifically, the authors demonstrate that the orbital time-evolution when expressed in the commonly used astrocentric elements can lead to spurious effects like a large scatter of the osculating semi-major axis and eccentricity of the outer planet about the secular value or a significant dependence of the integration results on the initial epoch. Unfortunately, the RV measurements of the HD 169830 system have not been published to the date. However, assuming that the best-fit parameters $`K,n,e,\omega ,T_\text{p}`$ (where, for every planet involved, $`K`$ is the semi-amplitude of RV variations, $`n`$ is the mean motion, $`e`$ is the eccentricity, $`\omega `$ is the periastron argument and $`T_\text{p}`$ is the time of periastron passage) describe a 2-Keplerian fit, they represent in fact the orbital parameters in the Jacobi coordinate system (Lee & Peale, 2003) and using them, we can recover the semi-major axes $`a`$ and the minimal masses of the companions. In some of the numerical experiments carried out in this paper, the inclination of the outer planet was changed by a tiny value, to 89.9, to allow the system enter the third spatial dimension for the purpose of the numerical integrations. The periastron passage of the inner planet is selected as the initial epoch of the osculating elements inferred from the 2-Keplerian solutions. The two ICs used in the computations are given in Table 1. The first, preliminary one, IC is labeled by ICI and the current one by ICII. In some experiments we modified ICI, by setting $`e_\text{c}10^6`$, to avoid numerical singularities. The elements of the HD 12661 system are given in Table 3 and we refer to them as IC0 throughout this work.

## 2 Numerical setup

We use the so-called fast indicators, MEGNO and FMA as the main numerical tools to study the orbital stability of the HD 169830 system. The indicator called MEGNO (the Mean Exponential Growth Factor of Nearby Orbits), has been invented by Cincotta & Simó (2000) and we have applied it in the studies of the exosystems’ dynamics described in a series of recent papers (for details and references on MEGNO see, for example, Goździewski et al., 2001; Goździewski & Maciejewski, 2003). The MEGNO is a very efficient tool—typically, it makes it possible to distinguish between regular and chaotic dynamics during an integration of the system carried out over $`10^4`$ orbital periods of the outermost planet. This indicator provides a direct estimate of the stability in the strict sense of the maximal Lyapunov Characteristic Number (LCN). As in the previous works, the MEGNO integrations are driven by a Bulirsh-Stoer integrator, namely the ODEX code (Hairer & Wanner, 1995). The relative and absolute accuracies of the integrator are set to $`10^{14}`$ and $`510^{16}`$, respectively. The zones of stability, corresponding to quasiperiodic motions of the system, are marked in the MEGNO maps by values close to 2. The positions of the nominal conditions are marked in the contour plots by the intersection of the two thin lines. In this work, we make yet another use of the integrations. During the MEGNO integrations we simultaneously analyze the osculating elements to study the short-term dynamics of the system. The osculating orbital elements feed the Frequency Map Analysis (FMA) algorithm (Laskar, 1993). The FMA, an already classic fast indicator, allows us to detect orbital resonances and determine diffusion rates of the fundamental frequencies in the planetary system (Robutel & Laskar, 2001). The diffusion rate is a quantity that cannot be determined by the LCN, however it also measures the lack of regularity of an orbit. Along a quasi-periodic solution the fundamental frequencies are constant. For a chaotic orbit these frequencies change and by calculating their diffusion rates, one can directly detect the macroscopic changes of the orbital elements. Thus by combining both methods, it is possible to derive extensive information on the system’s dynamics.

Since our computations concern mostly a short-term dynamics of the planetary system, the FMA is very useful in identifying the positions of the mean motion resonances (MMRs). In order to identify these resonances, during every integration of MEGNO, the complex functions $`f_k(t_i)=a_k\mathrm{exp}\text{i}\lambda _k(t_i)`$ are computed for planets $`k=\text{b,c}`$ at the discrete times $`t_i`$, where $`t_{i+1}t_i<0.5P_\text{b}`$, over the time $`2T`$ of the order of $`10^4P_\text{c}`$ ($`P_{\text{b,c}}`$ are the orbital periods of the inner and outer companion, respectively). Here, $`a_k`$ denotes the semi-major axis and $`\lambda _k`$ is the mean longitude of the planet. If the motion is quasiperiodic then in the time series $`\{f_k(t_i)\}`$, the frequency $`\nu _k`$ corresponding to the largest amplitude $`a_k^0`$ is one of the fundamental orbital frequencies of the system, called the proper mean motion, $`n_k`$ (Robutel & Laskar, 2001). Because the secular frequencies are much smaller than the orbital frequencies, in the first approximation the precession of the orbits can be neglected and MMRs can be identified through the condition $`q\nu _\text{b}p\nu _\text{c}0,`$ where $`p,q>0`$ are prime integers. The FMA gives the ratios $`\nu _\text{b}/\nu _\text{c}p/q`$, and then $`p`$ and $`q`$ can be easily resolved by the continued fraction algorithm.

After finding $`\nu _{\text{b,c}}^{(1)}`$ and $`\nu _{\text{b,c}}^{(2)}`$ where $`\nu _{\text{b,c}}^{(1,2)}`$ are the proper mean motions $`n_{\text{b,c}}`$ obtained after integrating the system over the time intervals $`[0,T]`$, and $`[T,2T]`$, respectively, the diffusion rate is determined through $`\sigma _k=1\nu _k^{(2)}/\nu _k^{(1)}`$ (Robutel & Laskar, 2001). Following these authors, if the motion is quasiperiodic, $`\sigma _k`$ is equal to zero or, computed with the FMA algorithm, is a very small quantity, typically less than $`10^8`$ while for an irregular chaotic motion, $`\sigma _k`$ is larger by several orders of magnitude. The inferred diffusion of the semi-major axis , can be derived through the Kepler law $`n^2a^3=\mu `$ as $`\mathrm{\Delta }n/n3/2(\mathrm{\Delta }a/a)`$. In this way, although the time span of the integrations is relatively very short, the most relevant instabilities of the motion can be detected and identified.

In this paper, we use the FMA code kindly provided by David Nesvorny on his web page333http://www.boulder.swri.edu/~davidn/fmft/fmft.html. According to the author, the code incorporates the idea of the Frequency Modified Fourier Transform, FMFT (Sidlichovsky & Nesvorny, 1997). We use a variant of the FMFT, with the so-called additional non-linear correction (see the cited paper for details). We have tested this code to make sure that it resolves the mean motions of non-interacting (Keplerian) orbits with the accuracy consistent with the internal precision of the method (set to $`10^{10}`$).

## 3 Orbital fits

The two distinctively different orbital solutions, ICI and ICII, announced within a short period of time by the Geneva group strongly suggest that at least one of them corresponds to a local minimum of $`(\chi _\nu ^2)^{1/2}`$. Unfortunately, it is not possible to verify their elements since the RV measurements are not available. In order to get at least some insight into the reliability of the newest set of orbital parameters, we have digitized the figure with RV measurements published on their WWW site 444http://obswww.unige.ch/~udry/planet/hd169830\_syst.html. The 93 data points we obtained differ from the real observations (for example, it is difficult to recover the exact moments of the observations), but they still properly describes the overall shape of the observed RV curve and its characteristic features. Note also that our set of ”measurements”555The set is available upon request from KG is smaller than the real one (112) because in some parts of the RV curve it was very difficult to clearly resolve all points due to a dense sampling of the real measurements. We assume that the observational errors are given by $`\sigma =\sigma _\text{0}+\sigma _\text{1}`$ where $`\sigma _\text{0}=8.5`$ m/s and $`\sigma _\text{1}`$ is an artificial normal random noise of 0 mean and the standard deviation about $`1`$, resulting in the ”errors” in the range $`[6,11]`$ m/s and roughly approximating the scale of the real errors.

We used the genetic algorithm scheme (GA) implemented by Charbonneau (1995) in his publicly available 666http://www.hao.ucar.edu/public/research/si/pikaia/pikaia.html code PIKAIA to to look for a global, coplanar orbital solution. The data were modeled with 2-Keplerian model (as in Goździewski et al., 2003) and the self-consistent Newtonian model, i.e., driven by the full $`N`$-body dynamics (Laughlin & Chambers, 2001). For the first model the GA code was executed a few hundreds of times and repeatedly found a solution very close to the best one, given in Table 2. Even though the number of our ”measurements” is smaller than in the real set, the similarity of the orbital solutions is very significant. Remarkably, the GA found very similar parameters for the outer planet (although the differences are larger than for the inner planet). The best fit parameters have the $`(\chi _\nu ^2)^{1/2}1.21`$ and the rms of $`10`$ m/s. This experiment was also repeated for two other data sets with 104 and 105 points (from two slightly different digitization) with the errors in the range $`[5,12]`$ m/s and $`[7,10]`$ m/s, respectively. The obtained best fit solutions were qualitatively the same as those given in Table 2 (let us note that the differences were most significant for $`e_\text{b}0.29`$, $`e_\text{c}0.29`$ and $`\omega _\text{c}270^{}`$). Our analysis indicates that the ICII by the Geneva group is really close to the global minimum of $`(\chi _\nu ^2)^{1/2}`$ and should not change at least if based on the currently available set of the RVs.

With the second model that accounts for the mutual interactions between planets, the best fit solution is not significantly better. It is characterized by the $`(\chi _\nu ^2)^{1/2}1.21`$ and the rms of $`9.5`$ m/s. The corresponding orbital parameters are very similar to those from the 2-Keplerian fit (see Table 2; for a reference, the synthetic RV-curves and the digitized ”data” points are also shown in Fig. 4). This result suggests that the current RV set does not allow to detect significant $`N`$-body coupling (for example, from low-order resonances). Finally, to roughly estimate the formal errors of the self-consistent Newtonian fit, we determined the confidence intervals of $`(\chi _\nu ^2)^{1/2}`$ (Press et al., 1992). The results are illustrated in Fig. 5. Clearly, the formal errors of $`a_\text{c}`$ and $`e_\text{c}`$ (semi-major axis and eccentricity of the outer planet) are substantial but the minimum of $`(\chi _\nu ^2)^{1/2}`$ is well localized. We note that to estimate the conservative errors in the best-fit parameters, the effect of stellar ”jitter” should be taken into account (as in Goździewski & Maciejewski, 2003).

## 4 Orbital stability of the HD 169830 system

Both orbital configurations of HD 169830 , ICI and ICII, appear to be stable. The evolution of the Jacobi osculating orbital elements computed over the time span of $`3.610^4`$ revolutions of the outer planet (about $`150,000`$ yr) for the ICI fit is illustrated in Fig. 1. During this time, the eccentricities (see Fig. 1c) vary with a large full-amplitude of about $`0.35`$. Fig. 1d reveals a presence of the secular apsidal resonance (SAR) in which the apsides are on the average antialigned. The semi-amplitude of the librations, $`\theta `$ (where $`\theta =\varpi _\text{b}\varpi _\text{c}`$ and $`\varpi _{\text{b,c}}`$ are the respective longitudes of periastron), is about $`90^{}`$. Apparently, the orbital elements vary in a regular way and this is strictly confirmed by the MEGNO signatures. The temporal changes of MEGNO, $`Y(t)`$, are shown in Fig. 1e and its mean value, $`Y(t)`$ (Fig. 1f), perfectly converges to 2. For a comparison, we computed these quantities for the HD 12661 system using the IC0 from Table 3. The results are shown in Figure 3. In both these cases, the variations of the orbital elements and the character of MEGNO convergence are qualitatively identical. However, the ICII fit for the HD 169830 system results in a qualitatively different orbital evolution (Fig. 2). The SAR does not seems to be present anymore. Nevertheless, the MEGNO signature (computed over about $`550,000`$ yr) indicates a stable, quasiperiodic configuration.

The MEGNO signatures have been verified by direct 1 Gyr integrations using the RMVS3 integrator from the SWIFT package Levison & Duncan (1994). We repeated the integrations with two different time steps, equal to 8 and 10 days. As one would expect, no instability occurs during this time and the orbital elements vary within the bounds determined by the short-term integrations. We do not show these results as they are basically an extension of the plots shown in Figures 1, 2 and 3. Obviously, due to the uncertainties of the orbital fits, such an examination of the isolated IC’s is not representative for the system’s dynamics. In order to find out whether the dynamics is robust to small adjustments of the initial condition, we have computed one-dimensional scans of $`Y`$ by changing the semi-major axis of the outer planet and keeping the other orbital parameters fixed at their initial values given in Table 1 for HD 169830 and 3 for HD 12661 . The results of this experiment are shown in Fig. 6. The MEGNO scan (the upper graphs of the respective panels in Fig. 6), computed with the resolution of about $`310^4`$ AU, reveals a number of spikes. Most of these spikes represent the MMRs between the planets. This identification is based on the FMA as described in the previous Section. We marked MMRs of the order $`p+q`$ not grater than about 20. These scans reveal that the dynamical environment of the nominal ICs is qualitatively the same for the ICI of the HD 169830 and the IC0 of HD 12661 systems. The outer planet lies in a stable zone between strong MMRs 6:1 and 7:1. The graph for the ICII of the HD 169830 system reveals a different dynamical environment. Clearly, the change in $`e_\text{c}`$ (compared to $`e_\text{c}`$ from ICI) ”shifted” the system to the zone between the 9:1 and 10:1 MMRs.

The FMA working with the MEGNO code enabled us to find out how sensitive both methods are to the instabilities of the motion. Bottom parts of the panels of Fig. 6 show the diffusion rate computed for the outer planet and estimated over its $`18,000`$ orbital periods. In the regions of quasiperiodic motions, as classified by MEGNO, the diffusion rate, $`\sigma _\text{c}`$, is smaller than about $`10^8`$. Such small values ensure us that the motion is close to a quasiperiodic evolution. In the MMRs zones $`\sigma _\text{c}`$ grows up to $`10^210^1`$. Hence in these areas, the osculating orbital elements exhibit macroscopic changes. Both algorithms are in an excellent accord, as they provide the same positions of the MMRS and very similar estimates of their widths. In other experiments, we noticed that for much shorter integration times, $`0.4510^4P_\text{c}`$, both algorithms can still detect all the relevant MMRs although in this case they do not allow to point out clearly enough some of the weak resonances. In fact, for the four times longer integration ($`1.810^4P_\text{c}`$), the FMA seems to be even more sensitive to the presence of weak resonances than the MEGNO algorithm is. We note here that by a compromise forced by the numerical efficiency of the code, in further runs of MEGNO, we have set the integration time to about $`0.910^4P_\text{c}`$.

## 5 Global dynamics of the HD 169830 system

To extend the above one-dimensional analysis, we can investigate the stability of the HD 169830 system in a few representative planes of its orbital parameters.

The program computing MEGNO simultaneously evaluated the maximal values of the eccentricities, $`e_{\text{b,c}}^{\text{max}}`$, attained during the integration time. We also stored the maximal value of semi-amplitude of the librations, $`\theta ^{\text{max}}`$, after every step (set to $`P_\text{c}`$) of renormalization of the variational equations. It helped us to detect the apsidal resonance and to estimate the semi-amplitude of the librations. The maximal value of the critical argument $`\theta `$ was taken relative to the center of the libration $`0^\mathrm{o}`$ or $`180^{}`$. To avoid the effects of a possible transition into the SAR, the determination of $`\theta ^{\text{max}}`$ was started after the first half of the integration period (about $`0.4510^4P_\text{c}`$). Finally, if $`\theta ^{\text{max}}<90^{}`$, then we treated this value as a semi-amplitude of the apsidal librations. The period of the integrations is relatively short but as we show in Goździewski (2003), such information on the short-term dynamics can still give us much insight into the global behaviour of the system. The results are illustrated in Figs. 79, where the left panels are for MEGNO, $`Y`$, the middle panels are for $`\theta ^{\text{max}}`$ and the right panels are for $`e_\text{b}^{\text{max}}`$. In these scans, the initial parameters that are varied are the coordinates of the maps and the other initial orbital elements are fixed at their nominal values given in Table 1.

Figure 7 is for the $`(a_\text{c},e_\text{c})`$-plane. Since the resolution of these maps is $`200\times 50`$ data points and the integration time is shorter than that of the one-dimensional $`a_\text{c}`$ scan (shown in Fig. 6), there is a lack of some fine resonance structures visible in Fig. 6. Yet the location of the dominant low-order MMRs: 6:1, 13:2 and 7:1 for ICI and 8:1, 9:1 and 10:1 for ICII, as well as their widths are clearly marked. The system would be chaotic for $`e_\text{c}`$ roughly grater than 0.5-0.6. In the central, stable parts of the MMRs, the maximal values of $`e_\text{b}`$ are small indicating their stabilizing influence on the motion. For ICI of the system around HD 169830 (as well as for IC0 of the HD 12661 system, not shown here), the plot for $`\theta ^{\text{max}}`$ reveals an extended zone of the SAR about the libration center of $`180^{}`$. For very small initial $`e_\text{c}0`$, the semi-amplitude of librations reaches the limiting $`90^{}`$. If the eccentricity becomes larger than $``$ 0.005 then the semi-amplitude decreases rapidly to $`60^{}`$-$`70^{}`$. The same type of $`(a_\text{c},e_\text{c})`$-scans of MEGNO for the HD 12661 system is shown in Goździewski (2003) (his Fig. 6) and Goździewski & Maciejewski (2003) (their Fig. 3). In both cases, the MEGNO structures around the nominal ICs are very similar. Interestingly, the $`(a_\text{c},e_\text{c})`$-maps of MEGNO for some of the fits to the RV data of the HD 160691 planetary system in Goździewski et al. (2003) (e.g., the fits GM5 and GM6, Fig. 5 in that paper), describe a dynamical setup very similar to these of the HD 169830 (with ICI) and HD 12661 systems. This is no longer true for the ICII of the HD 169830 system. In this case, the zone of the SAR is confined to a very small $`e_\text{c}`$ (but over wide range of $`a_\text{c}`$) and to the centers of the MMRs (see Fig. 7, bottom panels).

Figure 8 shows the behaviour of MEGNO, the semi-amplitude of apsidal libration and the maximal eccentricity of the inner planet for HD 169830 system on the $`(e_\text{b},e_\text{c})`$-plane. The nominal configuration of the system (denoted as previously by the intersection of the two thin lines) is located in an extended zone of quasiperiodic motions as seen in the MEGNO map (Fig. 8, panels in the left column). The $`e_\text{b}^{\text{max}}`$-maps (Fig. 8, right column) enable us to determine the border of the global instability established by these values of the initial $`e_\text{c}`$ which lead to large maximal eccentricities. For ICI, this border is substantially shifted towards larger initial eccentricities compared to the border of chaotic zone visible in the MEGNO map but overall it has a similar shape. The same effects are seen for ICII. In the SAR-maps (Fig. 8, middle column), for both IC’s one can find a very sharp border of this resonance in the region of small $`e_\text{b}`$. For ICI, it ends at the border of unstable motion present in the relevant $`e_\text{b}^{\text{max}}`$ panel. For ICII, the semi-amplitude of the librations is about $`60^{}`$$`70^{}`$ in this area but as for ICI (Fig. 7, middle panel), a very narrow strip of semi-amplitudes of $`90^{}`$ for the initial $`e_\text{c}0`$ is present. These large-amplitude librations for small $`e_\text{c}`$ are explained in the next section. Let us note that also in this test, the obtained picture of the stability zones for ICI of the HD 169830 system is qualitatively the same as for the HD 12661 system — see Fig. 5 and the discussion from Goździewski (2003). However, the SAR-map is quite different for ICII since the SAR is confined to a narrow zone of regular motion about small initial $`e_\text{c}`$ and also to unstable zone which is seen around $`(e_\text{b},e_\text{c})=(0.47,0.47)`$.

In the next set of experiments, we determine the zones of stability while the system inclination and hence planetary masses are varied. Assuming that the (unknown) inclinations of the orbits are the same $`i=i_\text{b}=i_\text{c}`$, the relative inclination $`i_{\text{rel}}`$ of the orbits depends on the difference of the nodal longitudes, $`\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Omega }_\text{c}\mathrm{\Omega }_\text{b}`$, according to the formula:

$$\mathrm{cos}(i_{\mathrm{rel}})=\mathrm{cos}(i_\text{b})\mathrm{cos}(i_\text{c})+\mathrm{sin}(i_\text{b})\mathrm{sin}(i_\text{c})\mathrm{cos}\mathrm{\Delta }\mathrm{\Omega }.$$

We assume that both orbits are prograde. By changing $`i`$, the masses are altered because the minimum masses, $`m_{\text{b,c}}\mathrm{sin}i`$, corresponding to edge-on orbits, must be preserved.

The results are shown in Fig. 9. For both ICI and ICII, even for large relative inclinations, the zones of stability extend up to very low system inclinations, $`i15^\mathrm{o}`$. This is supported by the $`e_\text{b}^{\text{max}}`$-scan which is shown in the right panel of this figure. In some unstable zones, seen in the MEGNO map for the ICI (Fig. 9, top-left panel), we can identify the 13:2 MMR. It quickly destabilizes the system in the areas about $`i45^{}`$—in these regions both eccentricities increases rapidly up to 0.8–0.9. In the strip about low inclinations, most likely thanks to the SAR, this effect is absent.

For ICI, due to the initial $`e_\text{c}0`$ and a geometrical singularity (for very small eccentricity the argument of periastron becomes undefined), our numerically computed $`\theta ^{\text{max}}`$-map (top-middle panel of Fig. 9) does not allow to clearly detect the zone of the SAR. However, in our preliminary study of the ICI we set $`e_\text{c}0.007`$ and then we found a wide zone of the SAR in the central part of the scanned area. The overall picture is qualitatively the same as in the case of the HD 12661 system (Goździewski & Maciejewski, 2003) where for low-$`i_{\text{rel}}`$ configurations, the SAR persists almost in the entire range of the system inclination $`i`$.

For ICII, the unstable zones in the MEGNO map are basically similar to those of the ICI of the HD 169830 and HD 12661 systems. Interestingly, all the maps for maximal $`e_\text{b}`$ reveal characteristic quarter-circle areas in which the eccentricity grows up to 1. In the strips near the central parts and corresponding to $`e_\text{b}1`$ the configuration becomes chaotic (compare the MEGNO- and $`e_\text{b}^{\text{max}}`$-maps). The analysis of the orbits in these areas shows that together with the excitation of the inner planet’s eccentricity, the inclination of this planet oscillates with a very large amplitude. It is illustrated in Figure 10 for the ICII and the initial relative inclination set to $`80^{}`$. Let us observe that in this case $`e_\text{b}`$ approaches 0.9 while the inclination goes up to $`140^{}`$. Nevertheless, the system still remains rigorously stable. This phenomenon seems to be an analogue of the Kozai resonance (Kozai, 1962) that is known for highly inclined asteroids and generally has been observed in triple systems in which the mass of the inner body is much smaller than the mass of the central and the outer body (Holman et al., 1997). In the multi-planetary exosystems, the masses of the perturber and the inner body are comparable which makes the problem more complex than in the asteroidal approximation. Let us note that a possibility of this phenomenon has been analyzed in $`\upsilon `$ Andromedae case (Chiang et al., 2001). The inclined configurations which we study here are rather special because it has been assumed that the initial inclinations are the same for both planets. However, a similar behaviour can be expected when the absolute inclinations are varied independently. It has been observed also in other hierarchical exosystems (see for example Goździewski, 2003, 2003). Since this phenomenon provides strong dynamical limits on the relative inclinations, it certainly deserves a detailed study which we plan to carry out in a future paper.

## 6 Secular apsidal resonance

The SAR in the HD 12661 exosystem was analyzed numerically in Goździewski (2003). We found an extended zone of the SAR in wide ranges of the semi-major axes, eccentricities and masses of the companions. Similar behavior of the SAR can be observed for the ICI of the HD 169830 system. It should be noted that those results were based on the IC of the HD 12661 system derived from the 2-Keplerian model by Fischer et al. (2003). Our attempts to analyze the same RV data set resulted in a rather different orbital solution Goździewski & Maciejewski (2003). Yet, the overall SAR features did not change. These results have been confirmed and greatly extended analytically by Lee & Peale (2003) in the framework of their octopole-level secular planetary theory.

Here, we literally quote that basics of their theory, which are relevant for our further analysis. This theory describes the secular dynamics in a coplanar hierarchical 2-planet system. It is assumed that MMRs are absent. The equations of motion originate from the Hamiltonian which is expanded up to the third order in parameter $`\alpha =a_\text{b}/a_\text{c}`$ and averaged over both mean longitudes $`l_j`$, $`j=\text{b,c}`$. The relevant parameters of the secular theory are defined through the momenta conjugated to the mean longitudes and arguments of periapse $`g_j=\omega _j`$:

$$\begin{array}{c}L_\text{b}=\frac{m_{}m_\text{b}}{m_{}+m_\text{b}}\sqrt{G(m_{}+m_\text{b})a_\text{b}},\hfill \\ L_\text{c}=\frac{(m_{}+m_\text{b})m_\text{c}}{m_{}+m_\text{b}+m_\text{c}}\sqrt{G(m_{}+m_\text{b}+m_\text{c})a_\text{c}},\hfill \\ G_j=L_j\sqrt{1e_j^2},\hfill \end{array}$$
(1)
where $`m_{}`$ is the mass of the parent star. The momenta are, respectively, the magnitude of the maximum possible angular momentum (for circular orbits) and the magnitude of the angular momentum. The sum of the remaining two momenta, $`H_j=G_j\mathrm{cos}i_j`$, is the $`z`$-component of the angular momentum of the system conjugate to the longitudes of the ascending node $`h_j=\mathrm{\Omega }_j`$. Let us note, that this formulation is valid in the general case of an inclined system. In coplanar configurations $`\omega _j`$ are undefined but then the longitudes of periastron $`\varpi _j`$ take over their role. In the averaged system, $`l_j`$ are absent in the Hamiltonian and $`L_j`$ are constants of motion. Because the octopole-level Hamiltonian depends only on the combination $`\varpi _\text{b}\varpi _\text{c}`$, it turns out that the total angular momentum $`G_\text{b}+G_\text{c}`$ is also the integral of motion in the averaged system. This reduces the coplanar dynamics to 1 degree of freedom, with $`e_\text{b}`$ or $`e_\text{c}`$, and $`\theta `$ as the relevant phase-space variables. The authors found that for the same constants:

$$\beta =\frac{5}{4}\frac{(m_{}m_\text{b})}{(m_{}+m_\text{b})}\alpha ,$$

$`\lambda =L_\text{b}/L_\text{c}`$, the same initial $`e_\text{b},e_\text{c}`$ and $`\theta =\varpi _\text{b}\varpi _\text{c}`$ the averaged equations of motion describe the same trajectories in the phase-space diagram of $`e_{\text{b,c}}(\theta )`$. For planetary masses much smaller than the central body, these trajectories and the amplitudes of eccentricity oscillations should be independent of the inclination of the system. Further, by introducing the parameter $`\gamma =(G_\text{b}+G_\text{c})/(L_\text{b}+L_\text{c})`$ which describes the non-dimensional total angular momentum of the system, one defines the critical value of $`\lambda `$,

$$\lambda _{\text{crit}}=\frac{2\gamma ^2}{53\gamma ^2}.$$

According to the equations of motion in the averaged system, if $`\lambda \lambda _{\text{crit}}`$ then libration of $`\theta `$ is almost certain with possibly large amplitude variations of both eccentricities. In general, the librations can take place about two libration centers, $`0^\mathrm{o}`$ or $`180^\mathrm{o}`$. Using the fits published in Fischer et al. (2003), the authors found that indeed, for HD 12661 system, $`\lambda \lambda _{\text{crit}}`$ and they identified an extended libration island of the SAR with the apsides antialigned in the exact resonance.

Following Lee & Peale (2003), we examined the ICs for the HD 169830 and HD 12661 systems by calculating these phase-space diagrams numerically. To obtain such diagrams, we fixed the constant level of the total angular momentum corresponding to the nominal IC and the phase curves were computed for varied $`e_\text{c}`$ ($`e_\text{b}`$ was calculated from the total momentum integral, the formula 1). The time of integrations was the same as for the derivation of the 2-dimensional stability maps. The results are shown in Figure 11. For both ICI of the HD 169830 and IC0 of the HD 12661 systems, the phase space is occupied by two large libration islands about $`\theta =0^{}`$ and $`\theta =180^{}`$ which is consistent with the secular theory. For the ICI of the HD 169830 system, $`\lambda =0.644`$ and $`\lambda _{\text{crit}}=0.9`$, so these parameters substantially differ. However, Fig. 11 shows that the occurrence of a libration mode still appears more likely than the $`\theta `$ rotations. For the IC0 of the HD 12661 , $`\lambda =0.75`$ is closer to its critical value $`\lambda _{\text{crit}}=0.87`$, and indeed, the the resonance island is more extended.

The configuration corresponding to ICI of the HD 169830 system is localized almost on the border of the island of the antialigned SAR. This border is a separatrix and if the initial $`e_\text{c}0`$, the semi-amplitude of $`\theta `$ librations approaches $`90^{}`$ as the system enters the vicinity of the separatrix. In the resonance island, the librations are about the libration center $`\theta =180^{}`$. The diagram also helps us to understand the large amplitudes of the eccentricities. Qualitatively, one obtains the same picture for the HD 12661 system. For the ICII of the HD 169830 system, we have $`\lambda =0.97`$ and $`\lambda _{\text{crit}}=0.77`$. In this case the nominal HD 169830 system lies in an extended area of $`\theta `$-rotations seen in Fig. 11.

The octupole theory allows us to quickly examine the extent of the SAR in the space of the orbital parameters. Because the elements of the inner planet are well determined, the overall picture of the HD 169830 dynamics depends mostly on the not too well constrained parameters of the outer planet. Assuming that strong mean motion resonances are absent, the secular dynamics can be investigated in detail through solving Equations (1)—(4) of Lee & Peale (2003). These equations describe the secular time-evolution of $`e_{\text{b,c}}`$ and $`\varpi _{\text{b,c}}`$. Unfortunately, there seems to be no explicit, analytical solution to these equations. Also a construction of a simple SAR criterion, like the one developed by Laughlin et al. (2002) or Ji et al. (2003) in the framework of the Laplace-Lagrange secular theory, does not appear to be possible, either. However, these equations are very simple for a numerical treatment and their direct integration is rapid, hundreds times faster than the $`N`$-body integrations. The SAR can be detected in the parameter space by looking for the maximal deviations of $`\theta `$ from the given libration center $`0^{}`$ or $`180^{}`$. It is the same method of detecting the SAR as the one applied in the full $`N`$-body integrations but here it is well justified because according to the generic secular dynamics, the SAR appears about two libration centers. Assuming that the integration time covers at least one secular period (so as it is not too short), we can easily detect the librations or rotations of the critical argument $`\theta `$.

Using this approach, we computed a number of the SAR-maps in the $`(a_\text{c},e_\text{c})`$ and $`(e_\text{b},e_\text{c})`$ domains for the initial conditions ICI, ICII and IC0. In order to directly compare the results, the ranges of the parameters were the same as for the numerical maps. The results are shown in Fig. 12. Generally, they agree in the regions corresponding to small eccentricities, for which both the shape of these regions and the semi-amplitude of librations are well reproduced by the secular theory. However, we can also find significant discrepancies between the analytical and numerical estimations. For example for ICII, in the $`(a_\text{c},e_\text{c})`$-plane one can find islands of the SAR which are absent in the map computed by solving the secular equations. The differences are also apparent in the $`(e_\text{b},e_\text{c})`$ plane for all examined IC’s. In these cases, the results agree only in the regions of relatively small $`e_\text{c}`$.

In order to understand these differences, we compared the $`(e_\text{c},\theta )`$ diagrams obtained by the direct numerical integrations and by solving the secular equations. The tests were performed for: (a) ICI in which we changed eccentricities to $`e_\text{b}=0.3,e_\text{b}=0.3`$, (b) ICII with $`a_\text{c}=3.51,e_\text{c}=0.47`$ (see bottom-middle panel in Fig. 7 corresponding to the center of the 9:1 MMR) (c) ICII with $`(e_\text{b}=0.47,e_\text{b}=0.47)`$ (corresponding to the island in the bottom middle panel in Fig. 8). The results are presented in Fig. 13. In all these cases, the full dynamics seems to be only roughly described by the secular approximation. The disagreement is obviously expected in the case (b) since the assumptions of the secular theory are violated (an MMR is present). In other two examples, the differences are most likely due to significant mutual interactions caused by large eccentricities resulting in deformations of the phase-space curves (case a) or a chaotic evolution (case c). Clearly the secular theory should be applied with caution.

## 7 Conclusions

In this paper we carry out a dynamical analysis of the recently discovered 2-planetary exosystem around HD 12661 . According to the preliminary initial orbital parameters announced by the Geneva Extrasolar Planet Search team (orbital solution ICI), this system resembled the other known planetary hierarchical exosystem HD 12661 and the HD 160691. The recently updated orbital solution (ICII) has been changed in a significant way and the close similarity of the HD 169830 and HD 12661 systems is not so apparent anymore.

The exosystems around HD 122661, HD 160691 and HD 169830 belong to a class of the hierarchical planetary systems. They are characterized by a relatively small ratio of the semi-major axes, $`0.10.3`$. Their dynamics should be analyzed in the Jacobi orbital elements (Lee & Peale, 2003; Goździewski et al., 2003). The commonly used astrocentric elements may lead to small shifts of the positions of the orbital resonances and to differences in the evolution of the osculating orbital elements. Following these papers we used the two orbital fits published by the discovery team but we recomputed the semi-major axes and masses of the companions. The current ICII was verified by our independent analysis of the approximate RV measurements which we obtained by digitizing a figure with the synthetic best-ft RV curve and the real observational points. Such an ”approach” to obtain the ”observations” was forced by the lack of access to the real data. Since the $`(\chi _\nu ^2)^{1/2}`$ function will typically have many local minima, we cannot be sure that a solution is the proper one without a global $`(\chi _\nu ^2)^{1/2}`$ analysis. As demonstrated on the solution ICI and ICII of the HD 169830 system, the difference between the global and a local minimum may lead to enormous differences in the overall dynamical picture of the system. Obviously, the digitization cannot provide very accurate moments of the observations, nevertheless such ”data” describe the the overall shape of the real RV measurements quite precisely and make it possible to obtain some insight into the quality of the best fit solution by the discovery team. Using the genetic algorithm to find the global minimum of $`(\chi _\nu ^2)^{1/2}`$ and describing the measurements with the 2-planet Keplerian and Newtonian models, we found almost the same orbital elements of the inner planet as these reported by the discovery team and very similar elements for the outer companion. The application of the 2-Keplerian and $`N`$-body models produced the same value of $`(\chi _\nu ^2)^{1/2}`$ and similar orbital solutions. It demonstrates the absence of strong interactions between the planets. It also favors the solution with a large separation between the companions rather than the previously announced configuration. Apparently, the first orbital solution was just a local minimum of the $`(\chi _\nu ^2)^{1/2}`$ function.

The two orbital solutions have different dynamical features. The ICI solution is close to an unstable 13:2 MMR and lies between 6:1 and 7:1 MMRs. The exosystem corresponding to ICII is located between 9:1 and 10:1 MMRs. Both configurations appear to be stable on the Gyr time scale as demonstrated by our long-term integrations, the fast indicator analysis and the secular octupole-level theory. In the phase space, both ICs are lie in the wide regions of stable motions in spite of large variations of the eccentricities (up to 0.5). Obviously, these conclusions may changed in the future as the time-span of the current RV measurements is shorter than the period of the outer planet. Moreover, the moderately constrained orbital elements of the outer planet do not allow us to exclude a proximity of the system to an unstable low-order resonance. The main dynamical difference between the two ICs is the SAR with the apsidal lines antialigned in the exact resonance. It is present for ICI in the wide ranges of the orbital parameters while for the ICII it seems to be ruled out. In the first case, the HD 169830 would be very similar to the hierarchical resonant systems. In this ”class” we can find the discussed HD 12661 and possibly the HD 160691 exosystem. The HD 37124 (Butler et al., 2003) system could also be a member of this group. However, in the HD 37124 system the apsides are most likely aligned (Ji et al., 2003). If the two outer planets are taken into consideration, the $`\upsilon `$ Andromedae (Butler et al., 1999) system could also be assigned to this group (Chiang & Murray, 2002). However, the new orbital solution ICII favors the proximity of the HD 169830 system to the HPS with the outer, very massive (possibly sub-stellar) companion like the HD 38529 (Fischer et al., 2003), HD 74156777http://obswww.unige.ch/~udry/planet/hd74156.html and HD 168443 (Marcy et al., 2001; Udry et al., 2002) systems. In these systems, the stability is maintained thanks to a relatively large separation of the companions. Nevertheless, due to very large masses involved and large amplitudes of the eccentricity variations, the mechanism providing the stability is puzzling.

The results of the numerical analysis of the orbital evolution and stability can be compared with the predictions by the secular octupole-level theory of Lee & Peale (2003). The results of both approaches are in accord for moderately low eccentricities of the outer planet. The discrepancies appear if both eccentricities are large in the regions of the phase space close to the MMRs and in the chaotic zones. It seems that although the secular theory makes it possible to explain the main features of the HD 169830 -like systems, the direct integrations are still necessary to understand the dynamics in detail.

For a better grasp on the dynamical picture of the HD 169830 system, we need to wait for some time, ideally until the observations cover a few orbital periods of the outer companion. Since the orbital elements of the inner companions are similar in the HD 169830 and HD 12661 systems and are already determined with high accuracy, there is still a chance that our preliminary hypothesis that the two planetary systems are dynamically related is plausible.

## 8 Acknowledgments

We thank Man Hoi Lee for a discussion about the Jacobi elements and Ji Janghui for useful remarks. This work is supported by the Polish Committee for Scientific Research, Grant No. 2P03D 001 22. M. K. is a Michelson Postdoctoral Fellow.


format: text

text:
# Resonant inelastic x-ray scattering study of holon-antiholon continuum in SrCuO2

## Abstract

We report a resonant inelastic x-ray scattering study of charge excitations in the quasi-one-dimensional Mott insulator SrCuO2. We observe a continuum of low-energy excitations, the onset of which exhibits a small dispersion of $`0.4`$ eV. Within this continuum, a highly dispersive feature with a large sinusoidal dispersion ($`1.1`$ eV) is observed. We have also measured the optical conductivity, and studied the dynamic response of the extended Hubbard model with realistic parameters, using a dynamical density-matrix renormalization group method. In contrast to earlier work, we do not find a long-lived exciton, but rather these results suggest that the excitation spectrum comprises a holon-antiholon continuum together with a broad resonance.

The separation of spin and charge degrees of freedom is one of the most important and fascinating properties of electrons in strongly correlated systems in one dimension. In particular, it is well known that in the one-dimensional (1D) Hubbard model, the low-energy physics is dominated by collective excitations of decoupled charge and spin degrees of freedom called holons and spinons, respectively Lieb68 . Experimentally, if one creates a hole by removing an electron, this hole is expected to decay into a spinon and a holon, which can be studied with angle resolved photoemission spectroscopy (ARPES) CKim . The situation is different for so-called “particle-hole” probes, such as optical spectroscopy, resonant inelastic x-ray scattering (RIXS), and electron energy loss spectroscopy (EELS). In these experiments, the total charge is conserved in the scattering process, so that an electron is simply moved from one site to another, creating a hole and a doubly occupied site. The decay of the hole creates a holon and a spinon, while the double-occupancy decays into an antiholon and a spinon. Since photons and electrons strongly couple to the charge sector, the behavior of holon-antiholon pairs can be studied with these particle-hole probes Neudert98 ; Hasan02 .

The so-called corner-sharing chain cuprates Sr2CuO3 and SrCuO2 are both charge-transfer insulators; that is, they have insulating gaps of $`2`$ eV arising from strong electron correlations. Since their crystal structure is highly anisotropic, the electronic structure remains 1D over a wide temperature range. Only at a very low temperature does magnetic order set in, due to the small interchain coupling ($`T_N2`$ K for SrCuO2 Matsuda97 ). Based on the commonly-used measure of quasi-one-dimensionality, $`T_N/J10^3`$, these compounds can be regarded as among the best realizations of quasi-1D systems Steiner76 . The quasi-1D nature of both SrCuO2 and Sr2CuO3 has been studied extensively with various experimental techniques, including magnetic susceptibility Motoyama96 , ARPES CKim , Raman scattering Misochko96 , optical spectroscopy Popovic01 , EELS Neudert98 , RIXS Hasan02 , and neutron scattering Zaliznyak03 .

In this Letter, we report a detailed study of the momentum dependence of the low-energy charge excitations in SrCuO2, utilizing the RIXS technique. We observe a continuum of excitations arising from the creation of particle and hole pairs. Within this continuum, a well-defined spectral feature with a large sinusoidal dispersion ($`1.1`$ eV) is observed. We have also measured the optical conductivity $`\sigma (\omega )`$ and carried out a dynamical density-matrix renormalization group (DDMRG) calculation of the extended Hubbard model. A comparison of the RIXS spectra with both the theoretical results and $`\sigma (\omega )`$ yields a consistent picture of charge excitations in this material: We find two dispersive spectral features, in contrast to the earlier RIXS study Hasan02 . The first feature is the onset energy of the holon-antiholon continuum, which exhibits small dispersion ($`0.4`$ eV). The second spectral feature is a more dispersive broad resonance that is the remnant of the strong-coupling exciton. This result contrasts earlier EELS work, in which the excitation spectrum was interpreted as being dominated by a bound exciton at large momenta Neudert98 .

In RIXS experiments, the incident x-ray energy is tuned to near the absorption edge of the particular element of interest, so that the inelastic scattering intensity of certain electronic excitations is resonantly enhanced. In addition, this gives an element specificity that is valuable in studying complex materials such as the cuprate compounds Hasan02 ; Hill98 ; Abbamonte99 ; Hasan00 ; LCO-PRL . The RIXS experiments were carried out at the Advanced Photon Source on the undulator beamline 9IDB with the same setup used in previous studies LCO-PRL . The scattering plane was vertical and the polarization of the incident x-ray was kept perpendicular to the Cu-O plaquette. Note that the corner-sharing CuO2 chain runs along the c-direction. A single crystal of SrCuO2 was grown using the traveling solvent floating zone method. The crystal was cleaved along the $`(\mathrm{0\hspace{0.33em}1\hspace{0.33em}0})`$ plane and mounted on an aluminum sample holder at room temperature in an evacuated chamber.

In order to find the resonance condition, we first carried out a study of the incident energy ($`E_i`$) dependence of the RIXS intensity. There are two features that show resonance behavior in the energy-loss ($`\omega E_iE_f`$) spectra, at $`\omega 3`$ eV and $`\omega 6`$ eV, respectively. The intensity of the 3 eV feature shows large enhancements at $`E_i8982`$ eV and $`E_i8990`$ eV, while the 6 eV feature is resonantly enhanced at $`E_i8990`$ eV and $`E_i8996`$ eV. We note that these three resonance energies correspond to the peaks in x-ray absorption spectra near the Cu K-edge Tranquada91 . The intensity at 6 eV likely arises from both the Cu K$`\beta _5`$ emission and a charge-transfer type excitation commonly found in the RIXS studies of cuprates Hill98 ; Abbamonte99 . Detailed studies of the $`E_i`$-dependence will be reported elsewhere, and here we fix $`E_i=8982`$ eV, and focus on the momentum dependence of the 3 eV feature.

Our main results are shown in Fig. 1, in which the RIXS intensity is plotted as a function of momentum and energy transfers. The pseudocolor scale of the observed intensity is shown on the right side in units of counts per second. The momentum transfer is shown in reduced units of $`q/2\pi `$ along the chain direction, and corresponds to the $`(\mathrm{0\hspace{0.33em}11}q)`$ position. Thus $`q/2\pi =0.5`$ is the 1D Brillouin zone boundary, while $`q/2\pi =0`$ is the zone center. The lower cut-off of $`\omega =0.5`$ eV shown in Fig. 1, is due to the large quasi-elastic contribution, and $`\omega =6`$ eV is the upper limit of accessible energy transfer in this setup. The plotted intensity is the raw RIXS intensity without any absorption corrections. Almost constant factors were obtained from three different normalization methods Norm , enabling us to compare the absolute intensities between different momentum transfers directly, except for $`q/2\pi =1`$. The anomalously low overall intensity at $`q/2\pi =1`$ arises from the very low angle of the incoming photons at this $`q`$.

The most prominent feature of the RIXS spectra shown in Fig. 1 is the highly dispersive feature at $`\omega 3`$ eV. The momentum dependence of the peak position of this feature, shown in Fig. 2, exhibits a clear sinusoidal dispersion with a bandwidth of 1.1 eV. However, closer inspection of Fig. 1 reveals the presence of additional spectral weight on the low-energy side around $`q/2\pi =0.5`$. This low-energy intensity around the zone boundary is evident for the $`q=\pi `$ data in Fig. 3, where the scattered intensity is plotted on a logarithmic scale. Since the momentum resolution at this position is about 0.1 reciprocal lattice units, it is unlikely that the additional low-energy intensity is due to a resolution effect. Thus, the best description of the observed spectra is that of a continuum of excitations, in which the dispersive feature resides. Note that this interpretation, particularly near $`q=\pi `$, is in contrast to previous EELS and RIXS results. In their EELS study of Sr2CuO3, Neudert et al. reported that the sharp spectral feature observed near $`q=\pi `$ is a bound exciton mode, due to a strong inter-site Coulomb interaction Neudert98 . The earlier RIXS study by Hasan et al. attributed the dispersive feature to the Mott gap edge; that is, a lower bound of particle-hole continuum Hasan02 . Our RIXS results suggest that the continuum starts at much lower energy, and that the particle-hole pairs do not form a bound exciton at $`q=\pi `$. The onset energy of the spectral weight extracted from each scan is shown in Fig. 3 as an arrow. The dispersion of these onset energies are plotted in Fig. 2. Note that the bandwidth of the onset energy dispersion ($`0.4`$ eV) is much smaller than that of the peak position, as discussed later.

In concert with the RIXS data, we have also studied the optical properties of SrCuO2, measuring the reflectivity of a cleaved surface between 6 meV and 6 eV, with the polarization along the c-axis. Optical conductivity, $`\sigma (\omega )`$ calculated from the reflectivity data using the Kramers-Kronig relation with an anchoring process, is plotted as a solid line in Fig. 4 and compared with the $`q=0`$ RIXS data. The validity of the Kramers-Kronig calculation was checked by directly measuring the complex dielectric function with an ellipsometry method. In the frequency range of interest ($`14`$ eV), identical spectra were obtained from both techniques. The $`q=0`$ RIXS data is expected to be different from that of $`\sigma (\omega )=\omega \mathrm{Im}[ϵ(0,\omega )]`$, since the response function measured in RIXS is believed to be close to a dielectric loss function, Im$`[1/ϵ(q,\omega )]`$. Instead, one should focus on the onset energy of the spectral weight, which should be independent of both the particular response function being measured and matrix-element effects. Both the RIXS data and the $`\sigma (\omega )`$ main feature exhibit onset energies around 1.7 eV, which corresponds to the continuum edge at the zone center.

To understand the observed continuum and the dispersive spectral feature, we have studied a simple one-band half-filled 1D extended Hubbard model (EHM) Stephan96 ; Tsutsui00 ; Essler01 ; Jeckelmann03 , in which the nearest-neighbor repulsion $`V`$ is explicitly included, in addition to the usual $`t`$ and $`U`$. It is clear that a full theoretical description of $`\mathrm{SrCuO}_2`$ requires a multi-band model based on the Cu $`3d`$ and O $`2p`$ orbitals. However, a precise and reliable determination of dynamical correlation functions in such models is still a formidable challenge. Several previous studies suggested that much of the essential physics of both $`\mathrm{Sr}_2\mathrm{CuO}_3`$ and $`\mathrm{SrCuO}_2`$ at low energies can be captured by the EHM Stephan96 ; Neudert98 ; Tsutsui00 ; Essler01 ; Jeckelmann03 . We have therefore used a DDMRG algorithm Jeckelmann03  to determine the optical conductivity and density response function in this model for realistic values of $`U`$ and $`V`$ on chains of up to 128 sites. All previous studies either dealt with unrealistic parameter values Neudert98  or suffered from severe finite-size effects caused by the very small system sizes amenable to exact diagonalization methods Hubsch01 .

The parameters $`t`$, $`U`$ and $`V`$ of the EHM were fixed by comparing DDMRG calculations of the optical conductivity $`\sigma (\omega )`$ to the experimental data and requiring compatibility with the value $`J0.23`$ eV for the Heisenberg exchange coupling $`J`$ obtained from neutron scattering Zaliznyak03 . As shown in Fig. 4, an excellent fit for the low-energy peak in $`\sigma (\omega )`$ is obtained for $`t=0.435`$ eV, $`V/t=1.3`$, $`U/t=7.8`$. For these values we obtain $`J0.24`$ eV from both $`1/U`$ expansions and DDMRG computations of the spinon dispersion in the single-particle spectral function. Note that the calculated $`\sigma (\omega )`$ in Fig. 4 agrees with the measured spectra on an absolute scale.

The excitation spectrum of the EHM with these parameters consists of scattering states of gapped, spinless charge $`e`$ collective modes – (anti)holons, and gapless, chargeless spin $`\pm 1/2`$ collective modes – spinons. In contrast to the strong-coupling results reported in Ref. Stephan96 , no excitonic excitations are found in any part of the Brillouin zone. The current operator is charge neutral and hence only excitations with equal numbers of holons and antiholons contribute to $`\sigma (\omega )`$. Most of the spectral weight of the low-energy peak in $`\sigma (\omega )`$ is attributed to excited states with one holon and one antiholon each. The contributions of the spin degrees of freedom to $`\sigma (\omega )`$ is small Essler01 .

The RIXS data is more difficult to model as at present there is no completely satisfactory theory of RIXS. However, it is clear that the RIXS process involves charge fluctuations that may be represented by the low energy excitations of the EHM Abbamonte99 ; Hasan02 ; Tsutsui00 . We have therefore determined the dynamical density-density correlation function $`𝒩(q,\omega )`$ by DDMRG. The results are shown in Fig. 3 as solid lines. At small momenta $`𝒩(q,\omega )`$ is proportional to $`q^2`$ and thus we use a normalization factor of $`\mathrm{sin}^2(q/2)`$ to compare $`𝒩(q,\omega )`$ at various momenta. We note that *on a lattice* the Fourier transform $`V(q)`$ of the Coulomb potential is approximately given by $`1/\mathrm{sin}^2(q/2)`$. Therefore, the renormalized $`𝒩(q,\omega )`$ in Fig. 3 is probably close to the dielectric loss function $`\mathrm{Im}[1/ϵ(q,\omega )]V(q)𝒩(q,\omega )`$ of the EHM. At $`q=0`$ there is a broad “band” of states with a peak around 2.3 eV. As we move towards the zone boundary, this feature narrows significantly. Importantly, however, even at $`q=\pi `$ the peak has an intrinsic width, that is, it is not a long-lived exciton mode. Furthermore, there is sizable intensity below the peak at $`q=\pi `$. This is in contrast to the strong-coupling $`U\mathrm{}`$ limit Stephan96 , where $`𝒩(q\pi ,\omega )`$ is dominated by an excitonic holon-antiholon bound state. Our interpretation is that as $`U`$ is reduced to smaller, more realistic values, the exciton acquires a finite lifetime and turns into a broad holon-antiholon resonance.

Let us try to relate these findings to the RIXS data. In Fig. 3, the normalized $`𝒩(q,\omega )`$ is directly compared with the RIXS spectra at equivalent $`q`$-positions. The onset of the RIXS spectrum shown in Fig. 3 seems to lie in the area where $`𝒩(q,\omega )`$ begins to show appreciable spectral weight. The onset energy of $`𝒩(q,\omega )`$ is plotted in Fig. 2 as a function of $`q`$ (crosses). The onset dispersion of $`𝒩(q,\omega )`$ follows the dashed line, which is the single-spinon dispersion relation shifted by a constant energy: $`\frac{\pi }{2}J|\mathrm{sin}(q)|+1.9`$. This spinon-like dispersion of the onset energy of $`𝒩(q,\omega )`$ is consistent with the low-energy field theory prediction Controzzi02 . Qualitative features of this spinon-like dispersion of the onset energy of $`𝒩(q,\omega )`$ appears to be consistent with features of the RIXS data (open circles, Fig. 2), such as the bandwidth of $`\pi J/2`$, and the $`q=\pi `$ point being a local minimum. In addition, we have compared the dispersion of the peak positions. As is indicated in Fig. 2, the peak dispersion in $`𝒩(q,\omega )`$ is very similar to that of the RIXS data. However, such a quantitative agreement between the RIXS peak positions and $`𝒩(q,\omega )`$ may well be coincidental, and further calculation of the full RIXS response function is necessary to draw any firm conclusions.

In summary, we have studied the energy and momentum dependence of low-energy charge excitations in quasi-1D SrCuO2, using both optical spectroscopy and resonant inelastic x-ray scattering. We observe a continuum of excitations which we associate with the creation of holon and antiholon pairs. The onset of the holon-antiholon continuum exhibits spinon-like dispersion. Within this continuum, a well-defined spectral feature with a large sinusoidal dispersion is also observed. By comparing the observed spectra with a DDMRG calculation of 1D EHM, we interpret this well-defined feature as a broad holon-antiholon resonance due to the inter-site Coulomb interaction. In contrast to earlier studies, we find that the holon-antiholon pairs do not form an exciton, and the continuum starts at much lower energy than previously reported.

We would like to thank F. Gebhard and C. C. Homes for invaluable discussions. The work at Brookhaven was supported by the US DOE, under contract No. DE-AC02-98CH10886. Use of the Advanced Photon Source was supported by the US DOE, under Contract No. W-31-109-Eng-38. H.B. acknowledges support by Optodynamics Center of the Philipps-Universität Marburg and thanks the Institute for Strongly Correlated and Complex Systems at BNL for hospitality and support.


text:
# Subthreshold Stochastic Resonance: Rectangular signals can cause anomalous large gains

## I Introduction

Most of the studies on the phenomenon of Stochastic Resonance (SR) in dynamical systems have been devoted to systems driven by sinusoidal terms (see PT96 ; RMP ; Wiesenfeld98 ; Anish99 ; Chemphyschem02  for reviews). Several analytical approximations have been put forward to explain SR. In the approach of McNamara and Wiesenfeld McNWie89 , the Langevin dynamics is replaced by a reduced two-state model that neglects the intra-well dynamics. The general ideas of Linear Response Theory (LRT) have been applied to situations where the input amplitude is small Hanggi78 ; PR82 ; RMP ; EPL89 ; GM89 ; PRA91 . In EPL89 ; PRA91  the Floquet theory has been applied to the corresponding Fokker-Planck description. For very low input frequencies, an adiabatic ansatz has been invoked PRA91 . Even though these alternative analytical approaches provide an explanation of SR for different regions of parameter values, their precise limits of validity remain to be determined. In recent work, EPLcasado ; FNLcasado  we have explored the validity of LRT for sinusoidal and multifrequency input signals with low frequency. Our results indicate a breakdown of the LRT description of the average behavior for low frequency, subthreshold amplitude inputs.

Several quantifiers have been used to characterize SR in noisy, continuous, systems. The average output amplitude, or the spectral amplification (SPA), has been studied in Refs. EPL89 ; PRA91  and the phase of the output average in Refs. dykman92 ; junhan93 , respectively. Those parameters as well as the signal-to-noise ratio (SNR) McNWie89 , exhibit a non-monotonic behavior with the noise strength which is representative of SR. An important quantity is the SR-gain, defined as the ratio of the SNR of the output over the input SNR. It has been repeatedly pointed out that the SR-gain can not exceed unity as long as the system operates in a regime described by LRT dewbia95 ; casgom03 . Beyond LRT there exists no physical reason that prevents the SR-gain to be larger than $`1`$, as it has been demonstrated in haninc00  for super-threshold sinusoidal input signals, and in analog experiments in ginmak00 ; ginmak01  for subthreshold input signals with many Fourier components and a small duty cycle casgom03 .

In this work, we will make use of numerical solutions of the Langevin equation following the methodology in casgom03 , to analyze SR in noisy bistable systems, driven by periodic piecewise constant signal with two amplitude values of opposite signs (rectangular signal) (see Fig. 1). There are several relevant time scales in the dynamics of these systems: i) $`t_{asym}`$, the time interval within each half period of the driving force, during which the diffusing particle sees an asymmetric constant two-well potential; ii) $`t_{inter}`$, the time scale associated with the inter-well transitions in both directions; and iii) $`t_{intra}`$, the time scale associated to intra-well dynamics. The inter-well and intra-well time scales depend basically on the noise strength $`D`$ and the amplitude of the driving term. The dependence of these two time scales with those parameters is certainly very different, being more pronounced for $`t_{inter}`$. Typically, for the range of parameter values associated with SR, the intra-well time scale is shorter than the inter-well one.

We will evaluate the long-time average behavior of the output and the second cumulant. These two quantities were studied some time ago by two of us Morillo95  for periodic rectangular driving signals. Here, we will further extend our work to the analysis of the correlation function and its coherent and incoherent parts. The knowledge of all these quantities provides a very useful information for the explanation of SR as indicated by the non-monotonic behavior with the noise strength of the output amplitude and the SNR. In particular, the knowledge of the incoherent part of the correlation function is of outmost importance for a correct determination of the SNR. Furthermore, for a given subthreshold amplitude, we will demonstrate that, if there exist a range of noise values such that $`t_{inter}`$ is shorter than $`t_{asym}`$, then it is possible to observe stochastic amplification and, simultaneously, SR-gains larger than unity. This is strictly forbidden by linear response theory, as we have recently shown casgom03 . Thus, a simultaneous appearance of stochastic amplification and SR-gains above $`1`$ implies a strong violation of linear response theory.

## II Model system and SR quantifiers

Let us consider a system characterized by a single degree of freedom, $`x`$, subject to the action of a zero average Gaussian white noise with $`\xi (t)\xi (s)=2D\delta (ts)`$ and driven by an external periodic signal $`F(t)`$ with period $`T`$. In the Langevin description, its dynamics is generated by the equation

$$\dot{x}(t)=U^{}\left[x(t)\right]+F(t)+\xi (t).$$
(1)
The corresponding linear Fokker-Planck equation (FPE) for the probability density $`P(x,t)`$ reads

$$\frac{}{t}P(x,t)=\widehat{}(t)P(x,t),$$
(2)
where

$$\widehat{}(t)=\frac{}{x}\left[U^{}(x)F(t)+D\frac{}{x}\right].$$
(3)
In the expressions above, $`U^{}(x)`$ represents the derivative of the potential $`U(x)`$. In this work, we will consider a bistable potential $`U(x)=x^2/2+x^4/4`$. The periodicity of the external driving $`F(t)`$ allows its Fourier series expansion in the harmonics of the fundamental frequency $`\mathrm{\Omega }=2\pi /T`$, i.e.,

$$F(t)=\underset{n=1}{\overset{\mathrm{}}{}}\left[f_n\mathrm{cos}(n\mathrm{\Omega }t)+g_n\mathrm{sin}(n\mathrm{\Omega }t)\right],$$
(4)
with the Fourier coefficients, $`f_n`$ and $`g_n`$, given by

$`f_n`$ $`=`$ $`{\displaystyle \frac{2}{T}}{\displaystyle _0^T}dtF(t)\mathrm{cos}(n\mathrm{\Omega }t),`$
$`g_n`$ $`=`$ $`{\displaystyle \frac{2}{T}}{\displaystyle _0^T}dtF(t)\mathrm{sin}(n\mathrm{\Omega }t).`$ (5)
Here, we are assuming that the cycle average of the external driving over its period equals zero. In this work, we will focus our attention to multi-frequency input forces with “rectangular” shape given by

$$F(t)=\{\begin{array}{cc}\hfill A:& 0t<\frac{T}{2}\hfill \\ \hfill A:& \frac{T}{2}t<T\hfill \end{array}$$
(6)
as sketched in Fig. 1. The external force remains constant at a value $`A`$ during each half-period and changes sign for the second half of the period. The duty cycle of the signal is defined by the time span the signal is nonzero over the total period of the signal; thus, the rectangular signal in Eq. (6) consequently possesses a duty cycle of unity.

The two-time correlation function $`x(t+\tau )x(t)_{\mathrm{}}`$ in the limit $`t\mathrm{}`$ is given by

$$x(t+\tau )x(t)_{\mathrm{}}=_{\mathrm{}}^{\mathrm{}}dx^{}x^{}P_{\mathrm{}}(x^{},t)_{\mathrm{}}^{\mathrm{}}dxxP_{1|1}(x,t+\tau |x^{},t),$$
(7)
where $`P_{\mathrm{}}(x,t)`$ is the time-periodic, asymptotic long time solution of the FPE and the quantity $`P_{1|1}(x,t+\tau |x^{},t)`$ denotes the two-time conditional probability density that the stochastic variable will have a value near $`x`$ at time $`t+\tau `$ if its value at time $`t`$ was exactly $`x^{}`$. It can been shown RMP ; PRA91  that, in the limit $`t\mathrm{}`$, the two-time correlation function $`x(t+\tau )x(t)_{\mathrm{}}`$ becomes a periodic function of $`t`$ with the period of the external driving. Then, we define the one-time correlation function, $`C(\tau )`$, as the average of the two-time correlation function over a period of the external driving, i.e.,

$$C(\tau )=\frac{1}{T}_0^Tdtx(t+\tau )x(t)_{\mathrm{}}.$$
(8)
The correlation function $`C(\tau )`$ can be written exactly as the sum of two contributions: a coherent part, $`C_{coh}(\tau )`$, which is periodic in $`\tau `$ with period $`T`$, and an incoherent part which decays to $`0`$ for large values of $`\tau `$. The coherent part $`C_{coh}(\tau )`$ is given by RMP ; PRA91

$$C_{coh}(\tau )=\frac{1}{T}_0^Tdtx(t+\tau )_{\mathrm{}}x(t)_{\mathrm{}},$$
(9)
where $`x(t)_{\mathrm{}}`$ is the average value evaluated with the asymptotic form of the probability density, $`P_{\mathrm{}}(x,t)`$.

According to McNamara and Wiesenfeld McNWie89 , the output SNR is defined in terms of the Fourier transform of the coherent and incoherent parts of $`C(\tau )`$. As the correlation function is even in time and we evaluate its time dependence for $`\tau 0`$, it is convenient to use its Fourier cosine transform, defined as

$$\stackrel{~}{C}(\omega )=\frac{2}{\pi }_0^{\mathrm{}}d\tau C(\tau )\mathrm{cos}(\omega \tau );C(\tau )=_0^{\mathrm{}}d\omega \stackrel{~}{C}(\omega )\mathrm{cos}(\omega \tau ).$$
(10)
The value of the output SNR is then obtained from:

$$R_{out}=\frac{lim_{ϵ0^+}_{\mathrm{\Omega }ϵ}^{\mathrm{\Omega }+ϵ}d\omega \stackrel{~}{C}(\omega )}{\stackrel{~}{C}_{incoh}(\mathrm{\Omega })}.$$
(11)
Note that this definition of the SNR differs by a factor $`2`$, stemming from the same contribution at $`\omega =\mathrm{\Omega }`$, from the definitions used in earlier works RMP ; PRA91 . The periodicity of the coherent part gives rise to delta peaks in the spectrum. Thus, the only contribution to the numerator in Eq. (11) stems from the coherent part of the correlation function. The evaluation of the SNR requires the knowledge of the Fourier components of $`C_{coh}(\tau )`$ and $`C_{incoh}(\tau )`$ at the fundamental frequency of the driving force. Thus, rather than knowledge of the entire Fourier spectrum, only two well defined numerical quadratures are needed. These are:

$$R_{out}=\frac{Q_u}{Q_l},$$
(12)
where

$$Q_u=\frac{2}{T}_0^Td\tau C_{coh}(\tau )\mathrm{cos}(\mathrm{\Omega }\tau ),$$
(13)
and

$$Q_l=\frac{2}{\pi }_0^{\mathrm{}}d\tau C_{incoh}(\tau )\mathrm{cos}(\mathrm{\Omega }\tau ).$$
(14)
The SNR for an input signal $`F(t)+\xi (t)`$ is given by

$$R_{inp}=\frac{\pi (f_1^2+g_1^2)}{4D}.$$
(15)
The SR-gain $`G`$ is consequently defined as the ratio of the SNR of the output over the SNR of the input, namely,

$$G=\frac{R_{out}}{R_{inp}}.$$
(16)
## III Numerical solution

Stochastic trajectories, $`x^{(j)}(t)`$, are generated by numerically integrating the Langevin equation \[Eq. (1)\] for every realization $`j`$ of the white noise $`\xi (t)`$, starting from a given initial condition $`x_0`$. The numerical solution is based on the algorithm developed by Greenside and Helfand hel79 ; grehel81  (consult also the Appendix in Ref. casgom03 ). After allowing for a relaxation transient stage, we start recording the time evolution of each random trajectory for many different trajectories. Then, we construct the long time average value,

$$x(t)_{\mathrm{}}=\frac{1}{N}\underset{j=1}{\overset{N}{}}x^{(j)}(t),$$
(17)
and the second cumulant,

$$x^2(t)_{\mathrm{}}x(t)_{\mathrm{}}^2=\frac{1}{N}\underset{j=1}{\overset{N}{}}\left[x^{(j)}(t)\right]^2\left[\frac{1}{N}\underset{j=1}{\overset{N}{}}x^{(j)}(t)\right]^2,$$
(18)
where $`N`$ is the number of stochastic trajectories considered. We also evaluate the two-time ($`t`$ and $`\tau `$) correlation function, i.e.,

$$x(t+\tau )x(t)_{\mathrm{}}=\frac{1}{N}\underset{j=1}{\overset{N}{}}x^{(j)}(t+\tau )x^{(j)}(t),$$
(19)
as well as the product of the averages

$$x(t+\tau )_{\mathrm{}}x(t)_{\mathrm{}}=\left[\frac{1}{N}\underset{j=1}{\overset{N}{}}x^{(j)}(t+\tau )\right]\left[\frac{1}{N}\underset{j=1}{\overset{N}{}}x^{(j)}(t)\right].$$
(20)
The correlation function $`C(\tau )`$ and its coherent part $`C_{coh}(\tau )`$ are then obtained using their definitions in Eqs. (8) and (9), performing the cycle average over one period of $`t`$. The difference between the values of $`C(\tau )`$ and $`C_{coh}(\tau )`$ allows us to obtain the values for $`C_{incoh}(\tau )`$. It is then straightforward to evaluate the Fourier component of $`C_{coh}(\tau )`$ and the Fourier transform of $`C_{incoh}(\tau )`$ at the driving frequency by numerical quadrature. With that information, the numerator and the denominator for the output SNR \[cf. Eqs. (12), (13) and (14)\], as well as the SR-gain \[cf. Eq. (16)\], are obtained.

## IV Results

### IV.1 Response to a rectangular driving force with fundamental frequency $`\mathrm{\Omega }=0.01`$

Consider an external driving of the type sketched in Fig. 1 with parameter values $`\mathrm{\Omega }=0.01`$, $`A=0.25`$. This amplitude is well below its threshold value defined, for each driving frequency, as the minimum amplitude that can induce repeated transitions between the minima of $`U(x)`$ in the absence of noise. For the input considered here, the threshold amplitude is $`A_T0.37`$. Note that this threshold value for the amplitude increases with increasing driving frequency.

In Fig. 2 we depict with several panels the behavior of the first two cumulants, $`x(t)_{\mathrm{}}`$ and $`x^2(t)_{\mathrm{}}x(t)_{\mathrm{}}^2`$, for several representative values of $`D`$ \[from top to bottom $`D=0.02`$ (panel a), $`D=0.04`$ (panel b), $`D=0.06`$ (panel c), $`D=0.1`$ (panel d) and $`D=0.2`$ (panel e)\]. Notice that due to the transients, the time at which we start recording data, $`t=0`$ in the graphs, does not necessarily coincide with the start of an external cycle. The average is periodic with the period of the driving force, while the second cumulant, due to the reflection symmetry of the potential EPLcasado , is periodic with a period half of the period of the forcing term.

Next we consider the case of small noise intensity $`D`$ (say, $`D=0.02`$ as in panel a). The noise induces jumps between the wells. In each random trajectory, a jump between the wells has a very short duration, but the instants of time at which they take place for the different stochastic trajectories are randomly distributed during a half-cycle. At this small noise strength, the jumps are basically towards the lowest minimum. Thus, because of this statistical effect, the average behavior depicts the smooth evolution depicted in panel a of Fig. 2, without sudden transitions between the wells. The evolution of the second cumulant adds relevant information. The fact that it is rather large during most of a period indicates that the probability density $`P_{\mathrm{}}(x,t)`$ is basically bimodal during most of the external cycle. It is only during very short time intervals around each half-period that the probability distribution becomes monomodal around one of the minima, and, consequently, the second cumulant is small. The bimodal character arises from the fact that the noise is so small in comparison with the barrier heights that jumps over the barrier are rather infrequent during each half a period.

As the noise strength increases, the time evolution of $`x(t)_{\mathrm{}}`$ follows closely the shape of the external force, (cf. see panels b, c, d in Fig. 2 for $`0.04D0.1`$). This behavior indicates that, for these parameter values, the jumps in the different random trajectories are concentrated within short time intervals around the instants of time at which the driving force switches sign. The second cumulant remains very small during most of a period, except for short time intervals around the switching times of the external driver. Thus, for these intermediate values of $`D`$, the probability distribution, $`P_{\mathrm{}}(x,t)`$ is basically monomodal, except for small time intervals around the switching instants of time of the periodic driver. Finally, as the noise strength is further increased ($`D>0.1`$), the probability distribution remains very broad most of the time. Even though a large majority of random trajectories will still jump over the barrier in synchrony with the switching times of the driver, the noise is so large in comparison with the barrier heights, that the probability of crossings over the barrier in both directions can not be neglected at any time during each half cycle. The probability distribution remains bimodal during a whole period, but asymmetric: the larger fraction of the probability accumulates around the corresponding lower minima of the potential during each cycle. Therefore, the average output amplitude decreases, while the second cumulant depicts plateaus at higher values than for smaller noise strengths (compare panels e and c in Fig. 2).

In Fig. 3, we plot the coherent (left panels), $`C_{coh}(\tau )`$, and incoherent (right panels), $`C_{incoh}(\tau )`$, components of the correlation function $`C(\tau )`$ for the same parameter values as in Fig. 2. The coherent part shows oscillations with a period equal to that of the driving force. Its shape changes with $`D`$. The amplitude of the coherent part does not grow monotonically with $`D`$. Rather, it maximizes at $`D0.06`$, which is consistent with the observed behavior of $`x(t)_{\mathrm{}}`$ in Fig. 2. This is expected as the evaluation of $`C_{coh}(\tau )`$ involves only the time behavior of $`x(t)_{\mathrm{}}`$ at two different instants of time.

Two features of the behavior of $`C_{incoh}(\tau )`$ are relevant: its initial value and its decay time. The initial value of the incoherent contribution, $`C_{incoh}(0)`$, is given by the cycle average of the second cumulant. It has a non-monotonic behavior with $`D`$. For $`D=0.02`$, $`C_{incoh}(0)`$ is large, consistent with the fact that the second cumulant at this noise strength is appreciably different from $`0`$ during a substantial part of a period. As the value of $`D`$ increases, ($`D0.1`$), the value of $`C_{incoh}(0)`$ decreases. This is expected as the second cumulant is large just during those small time intervals where most of the forward transitions take place every half-period. For still larger values of $`D`$ there are frequent forward and backward jumps that keep the stationary probability bimodal, and therefore, the initial value $`C_{incoh}(0)`$ increases. For $`D=0.02`$, the decay of $`C_{incoh}(\tau )`$ is very slow, although the decay time is still shorter than the duration of half a period of the driving force. As $`D`$ increases, the decay time of the incoherent part becomes shorter. It is worth to point out that the intra-well noisy dynamics manifests itself in the behavior of $`C_{incoh}(\tau )`$. This is most clearly confirmed noticing the fast initial decay observed in panel d. For smaller values of $`D`$, this feature is masked by the long total relaxation time scale, while for $`D=0.2`$, the noise strength is so large that there is not a clear-cut separation between inter and intra-well time scales.

The above considerations allow us to rationalize the behavior of the several quantifiers used to characterize SR. Their behaviors with $`D`$ for $`A=0.25`$, $`\mathrm{\Omega }=0.01`$ are depicted in Fig. 4. It should be noticed that the lowest value of the noise strength used in the numerical solution of the Langevin equation is $`D=0.01`$. For this noise strength, the values of $`Q_u`$ and $`R_{out}`$ are very small, although not zero. For even lower noise strengths the task becomes computationally very demanding and expensive, due to the extremely slow decay of the correlations. For $`D`$ sufficiently small, however, one does expect $`Q_u`$ to be larger than $`Q_l`$, and, consequently, an increase of the numerical $`R_{out}`$ as $`D`$ is lowered.

The quantity $`Q_u`$ defined in Eq. (13) depicts a non-monotonic dependence on $`D`$ typical of the SR phenomenon. Its behavior is expected from the dependence of the amplitude of $`C_{coh}(\tau )`$ with $`D`$ in Fig. 3.

A non-monotonic behavior with $`D`$ for the numerically evaluated $`Q_l`$ is also observed. The initial value, $`C_{incoh}(0)`$ and the decay time of $`C_{incoh}(\tau )`$ are important in the evaluation of $`Q_l`$ \[see Eq. (14)\]. For $`D=0.01`$, the decay time of the incoherent part of the correlation function is longer than half a period of the driving force, while for $`D=0.02`$, it is somewhat shorter than $`T/2`$. Consistently with Eq. (14), the value of the integral for $`D=0.01`$ is smaller than for $`D=0.02`$. As $`D`$ is further increased, the influence of the cosine factor in Eq. (14) becomes less important as the decay time is much shorter than $`T/2`$. The drastic fall in the $`Q_l`$ values observed for $`0.02<D<0.1`$ is due to the decrease of $`C_{incoh}(0)`$ with $`D`$ (see panels a-c in Fig. 3) and the shortening of the decay time. As $`D`$ is increased further, $`C_{incoh}(0)`$ increases and, consequently, $`Q_l`$ also increases slightly.

Taking into account the definition of the SNR, \[cf. Eq. (12)\], its behavior with $`D`$ is not surprising. The numerically obtained SNR peaks at $`D=0.08`$, a slightly different value of $`D`$ from the one at which $`Q_u`$ peaks.

### IV.2 Anomalous SR-gain behavior for subthreshold driving

The SR-gain is defined in Eq. (16). The numerically determined SR-gain shows a most interesting feature: we observe a non-monotonic behavior versus $`D`$, with values for the gain exceeding unity (!) for a whole range of noise strengths. This is strictly forbidden within LRT casgom03 ; therefore, the fact that the SR-gain can assume values larger than unity reflects a manifestation of the inadequacy of LRT to describe the system dynamics for the parameter values considered.

To rationalize this anomalous SR-gain behavior, we notice that the role of the noise in the dynamics is twofold. On the one hand, it controls the decay time of $`C_{incoh}(\tau )`$. On the other hand, the noise value is relevant to ellucidate whether the one-time probability distribution is basically monomodal or bimodal during most of the cycle and, consequently, it controls the initial value $`C_{incoh}(0)`$. As discussed above, if $`D`$ is small, the decay time is very large compared to $`t_{asym}`$, and the one-time probability distribution is essentially bimodal. For large values of $`D`$, the decay of $`C_{incoh}(\tau )`$ is fast enough, and the distribution is also bimodal. The large SR-gain obtained here requires the existence of a range of intermediate noise values such that: i) $`C_{incoh}(\tau )`$ decays in a much shorter time scale than $`t_{asym}`$ and, ii) the one-time probability distribution remains monomodal during most of the external cycle.

### IV.3 Response to a rectangular input driver with fundamental frequency $`\mathrm{\Omega }=0.1`$

As mentioned before, there are several time scales which are important for the phenomenon of stochastic amplification and gain. In the previous subsection we have considered an input frequency small enough so that the inequality $`t_{asym}>t_{inter}`$ holds for a range of noise values. Next, we shall analyze the system response to a driving force with fundamental frequency $`\mathrm{\Omega }=0.1`$, ten times larger than in the previous case. We will take the same input amplitude as in the previous case, $`A=0.25`$, which is still subthreshold. For this input frequency, the threshold value for the amplitude strength is determined by numerically solving the deterministic equation to yield $`A_T0.39`$.

The behavior of the first two cumulants for several values of the noise strength is depicted in Fig. 5 \[from top to bottom $`D=0.02`$ (panel a), $`D=0.04`$ (panel b), $`D=0.06`$ (panel c), $`D=0.1`$ (panel d) and $`D=0.2`$ (panel e)\]. For all values of $`D`$, the second cumulant remains large for most of each half-period. By contrast with the lower frequency case, we detect no values of $`D`$ for which the probability distribution is monomodal for a significant fraction of each half period.

In Fig. 6 the behavior of the coherent (left panels) and incoherent (right panels) parts of the correlation function are presented for (from top to bottom) $`D=0.02`$ (panel a), $`D=0.04`$ (panel b), $`D=0.06`$ (panel c), $`D=0.1`$ (panel d) and $`D=0.2`$ (panel e). The amplitude of the coherent oscillations shows a nonmonotonic behavior with $`D`$. The incoherent part has initial values which remain very large in comparison with the corresponding ones for $`\mathrm{\Omega }=0.01`$ (compare with Fig. 3), consistently with the large value of the second cumulant. The decay times are roughly the same for both frequencies.

In Fig. 7 we show the behavior of the several SR quantifiers as a function of $`D`$. The comparison of Figs. 4 and 7 indicates that $`Q_u`$ , $`Q_l`$ and $`R_{out}`$ have the same qualitative behavior for both frequencies. The non-monotonic dependence on $`D`$ of $`Q_u`$ and $`R_{out}`$ are indicative of the existence of SR (for both frequencies) for the subthreshold input amplitude and in the ranges of $`D`$ values considered. The most relevant quantitative difference is that for $`\mathrm{\Omega }=0.1`$ the SR-gain remains less than unity.

## V Conclusions

With this work, we have analyzed the phenomenon of SR within the context of a noisy, bistable symmetric system driven by time periodic, rectangular forcing possessing a duty cycle of unity. The numerical solution of the Langevin equation allows us to analyze the long-time behavior of the average, the second cumulant and the coherent and incoherent parts of the correlation function. For subthreshold input signals we determined the SNR, together with its numerator and denominator evaluated separately, for a wide range of noise strengths $`D`$.

As a main result we find the simultaneous existence of a typical non-monotonic behavior versus the noise strength $`D`$ of several quantifiers associated to SR; in particular SR-gains larger than unity are possible for a subthreshold rectangular forcing possessing a duty cycle of unity. This finding is at variance with the recent claim in Refs. ginmak01 ; makging02  that pulse-like signals with small duty cycles are needed to obtain SR-gains larger than unity. This unexpected result occurs indeed whenever the inequality $`t_{asym}>t_{inter}`$ holds. This is most easily achieved with low frequency inputs. As the input frequency increases, that inequality is not satisfied for sufficiently small values of $`D`$: even though SR then still exists, it is not accompanied by SR-gains exceeding unity. Furthermore, in Refs. ginmak01 ; makging02 , SR-gains larger than unity are only obtained with input amplitudes larger than $`0.8A_T`$. By contrast, in this work, we have shown that such a large value for the input amplitude is not needed (we have used $`A0.68A_T`$).

The simultaneous occurrance of SR and SR-gains larger than unity is associated to the fact that, for some range of noise values, the decay time of the incoherent part of the correlation function is much shorter than $`t_{asym}`$ and also the probability distribution is basically monomodal during most of the cycle of the driving force.

###### Acknowledgements.

We acknowledge the support of the Dirección General de Enseñanza Superior of Spain (BFM2002-03822), the Junta de Andalucía, the DAAD program ”Acciones Integradas” (P.H., M. M.) and the Sonderforschungsbereich 486 (project A10) of the Deutsche Forschungsgemeinschaft.