text:
# 1 Introduction

## 1 Introduction

Because pulsars emit the majority of their spin-down energy in gamma-rays, understanding their emission at these energies is crucial for forming a complete picture of pulsar energetics. However, while the radio pulsar population now numbers almost 1500, pulsed gamma rays have been detected from less than 10 sources. With these sparse statistics, addressing important issues such as how gamma-ray luminosity depends on spin-down parameters, the gamma-ray pulsar emission mechanism, the relationship between radio and gamma-ray beams, the pulsar contribution to the unidentified EGRET source population, and the pulsar science prospects of future gamma-ray missions is difficult. For this reason, we developed a likelihood analysis which uses the EGRET pulsar detections, upper limits, and diffuse background measurements to characterize some properties of the gamma-ray pulsar population. In this paper, we outline our updates to the analysis of McLaughlin & Cordes 2000 (hereafter MC00), present the new results, and discuss the prospects of AGILE and GLAST for pulsar science. We also discuss issues involved in searching for gamma-ray pulsars and present the methodology and results from searches for gamma-rays from new radio pulsars. Such searches are difficult due to the sparseness of gamma-ray photons and the lack of contemporaneous radio ephemerides.

## 2 Model

Our likelihood function is a product of the individual likelihoods for pulsar detections, upper limits and diffuse background measurements. We model a pulsar’s gamma-ray luminosity $`L`$ as a power-law in period $`P`$ (in seconds) and magnetic field $`B_{12}`$ (in units of $`10^{12}`$ G) such that $`L=\gamma P^\alpha B_{12}^{}{}_{}{}^{\beta }`$. We assume a spin-down law with braking index $`n`$ to calculate a population-averaged luminosity. To calculate all likelihoods we assume a broad beaming solid angle of $`2\pi `$. To calculate the diffuse background likelihood, we assume a constant pulsar birthrate of 1/100 yr, a Galactic age of $`10^{10}`$ years, and a maximum gamma-ray efficiency of 1/2. We allow pulsars to contribute up to 10% of the total diffuse flux and assume that they are distributed in a Gaussian disk of scale 6 kpc with exponential halo of scale 0.5 kpc and molecular ring at 4 kpc with width 1.5 kpc. Please see MC00 for detailed descriptions of the likelihood calculations.

In MC00, we used the Taylor & Cordes (1993) electron density model to calculate pulsar distances. In this analysis, we use Cordes & Lazio (2002). This has not changed the results dramatically, but has been important for some individual objects. For example, the implied gamma-ray efficiency for B1055$``$52 has decreased from 30% to 5%. In MC00, we assumed a single value for the population’s magnetic field $`B`$ and initial spin period $`P_0`$. In this analysis, we fit for a lognormal distribution for $`B`$ characterized by mean $`<B>`$ and rms $`\sigma _B`$ and a flat distribution for $`P_0`$ bounded by $`P_{0,min}`$ and $`P_{0,max}`$.

## 3 Results

Our EGRET data (described in MC00) consist of 7 pulsar detections, 353 pulsed flux upper limits and 3 diffuse background measurements. We calculate the likelihood for a range of values for $`\alpha `$, $`\beta `$, $`\gamma `$, $`n`$, $`<B>`$, $`\sigma _B`$, $`P_{0,min}`$, and $`P_{0,max}`$ and find well-defined maxima in $`\alpha `$, $`\beta `$, and $`\gamma `$ of 1.7, 1.2, and 32.4. We cannot constrain $`n`$, $`<B>`$, $`\sigma _B`$, $`P_{0,min}`$, or $`P_{0,max}`$ as for any reasonable values the pulsar contribution to the diffuse flux is $``$ 10%. As the best-fit luminosity law therefore only depends on detections and upper limits, it is independent of all assumptions aside from beaming solid angle and distance. Table I shows that our best-fit law is most similar to ones in which the luminosity is proportional to the voltage drop across the polar cap. It is not consistent with luminosity proportional to spin-down energy and is also quite different from the best-fit OSSE (McLaughlin et al. 2000) and radio (Arzoumanian et al. 2002) laws.

We can use the best-fit law and the assumptions of Section 2 to calculate the gamma-ray flux distribution for a model population of pulsars (see Figure 1). To do this, we must assume values for $`n`$ and for the $`B`$ and $`P_0`$ distributions. For only five pulsars has $`n`$ been measured, with values ranging from 1.5 to 3 (see Zhang et al. 2001 and references therein). We adopt $`n=2.5`$, as measured reliably for the Crab (Lyne et al. 1993). Fitting a lognormal distribution to the magnetic fields (i.e. $`B_{12}^{}{}_{}{}^{2}=10^{12}P\dot{P}`$) of 1244 non-recycled pulsars, we find a mean of 12.1 G and rms of 0.65 G. $`P_0`$ has only be estimated for eight pulsars associated with supernova remnants with independent age estimates and/or measurable pulsar proper motions. With estimates ranging from $`<14`$ ms for J0537$``$6910 to 139 ms for J0538+2817 (see Migliazzo et al. 2002, Kramer et al. 2003, and references therein), we adopt a flat distribution from 10 to 150 ms.

This analysis shows that, given EGRET’s point source sensitivity of $`10^7`$ ph cm$`^2`$s-1 $`10^{32.9}`$ ergs s-1 kpc-2 (for a spectral index of $`2`$), EGRET should have detected roughly 40 pulsars as point sources. For a pulsed-search sensitivity of $`6\times 10^7`$ ph cm-2s-1 $`10^{33.7}`$ ergs s-1 kpc-2 (for a 2-week integration, a duty cycle of 1/2, and a minimum $`N_s^2/N_t`$ of 50, where $`N_s`$ is the number of source counts and $`N_t`$ is the number of total counts), EGRET could have detected six pulsars in blind searches. Of the 170 unidentified EGRET sources, 73 are within 10 degrees of the Galactic plane and 47 are non-variable (McLaughlin 2001), suggestive of a pulsar origin and consistent with our estimate of 40 pulsars detected as point sources. We can also compare with the analysis of Kramer et al. (2003) which found that, of the 48 positional coincidences between radio pulsars and unidentified EGRET sources, $`19\pm 6`$ are likely to be real. The estimate of 40 EGRET pulsars is dependent upon a beaming solid angle of $`2\pi `$. For a beaming solid angle of $`\pi `$, EGRET should have detected 20 pulsars as point sources. However, given the analysis of Kramer et al., this would imply that all pulsars detected by EGRET are currently known radio pulsars, a hypothesis which seems unlikely given that we know of at least one gamma-ray pulsar (Geminga) which is not detected at radio wavelengths.

In Table II, we list all radio pulsars with high model-predicted fluxes that are associated with unidentified EGRET sources111J1747–2958 actually lies outside of the 95% position confidence contour of 3EG J1746–2851. However, because of diffuse model and positional uncertainties, especially towards the Galactic center, we believe this is a possible association. . Six of these 12 pulsars have been detected in the Parkes Multibeam Pulsar Survey, which has discovered half of the 80 radio pulsars with ages $`<`$ 100 kyr. We also list the EGRET fluxes (in units of $`10^{33}`$ ergs s-1 kpc-2) and variability indices (McLaughlin 2001) of these sources. While all of the known EGRET pulsars have non-variable fluxes, some of these unidentified EGRET sources appear to be highly variable. It is difficult to determine whether these are spurious coincidences or if there is indeed some method through which pulsars may produce variable gamma-ray emission. One method is through searching for gamma-ray pulsations (see Section 4).

Given the expected sensitivities of AGILE and GLAST, we can predict the pulsar populations that each instrument will see. In Table III, we list the sensitivities (in units of $`10^{31}`$ ergs s-1 kpc-2), and the numbers of expected detections of pulsars as point sources and in blind periodicity searches. In column two, the number of known radio pulsars that are expected to be detected is listed in parentheses (i.e. for EGRET, our model predicts that 17 out of the 40 total pulsar detections are known radio pulsars). In column four, the number in parentheses lists how many EGRET unidentified sources could be identified as pulsars in a blind periodicity search. These are optimistic estimates as they assume a 2-week, on-axis exposure with fresh gas. Note that GLAST will be able to detect pulsations from ALL of the unidentified sources in blind periodicity searches, unambiguously determining their nature. All of the estimates in Table III depend on a beaming solid angle of $`2\pi `$ and will decrease linearly with decreasing solid angle.

## 4 Gamma-Ray Pulsar Searching

With many pulsars expected to be detected at gamma-ray energies but not at radio wavelengths, developing methods to search for gamma-ray periodicities is crucial. We are applying such methods to search EGRET data for gamma-rays from the new radio pulsars listed in Table II. While we have current radio ephemerides for these objects, they are too imprecise to allow extrapolation back to the time the data were taken. Therefore, searching over a wide range of frequency $`f`$ and frequency derivative $`\dot{f}`$, centered around those predicted by the current ephemeris, is necessary. To calculate the range to search, we assume errors in $`f`$ and $`\dot{f}`$ of roughly $`10^8`$ s-1 and $`10^{12}`$ s-2, typically many times greater than the quoted timing errors, to account for glitches and timing noise. We space trial values of $`f`$ and $`\dot{f}`$ so that the smearing in the final folded profile is $`<10`$% of the period. For each trial, we barycenter all photons to the pulsar’s position and fold all those with energies $`E_\gamma 100`$ MeV within an energy-dependent acceptance cone of radius $`\theta 5^o.85(E_\gamma /100\mathrm{M}\mathrm{e}\mathrm{V})^{0.534}`$, weighting photons by the PSF (e.g. Ramanamurthy et al. 1996). We compute the significance of each profile using $`\chi ^2`$, $`H`$ (e.g. de Jager 1994), and $`Z^2`$ (e.g. Buccheri et al. 1983), normalizing by the number of trials, and inspect the most significant signals as determined by all of those three methods.

We have applied this method to all of the pulsars listed in Table II and to the known EGRET pulsars B0531+21 and B1706–44. For these known pulsars, we use ephemerides based on only 2 yrs of radio data for comparison with the timing solutions of similar timescale for the new pulsars. Examining the known pulsar results show that our inclusion of all photons above 100 MeV is optimal for detection; a lower cutoff results in less significant detections due to increased background photons but higher cutoffs result in too few pulsar photons. We find that for both known pulsars and in all viewing periods (VPs), weighting by the PSF increases the significance of detections. We find that $`\chi ^2`$, $`H`$, and $`Z^2`$ generally return the same values for best $`f`$ and $`\dot{f}`$ but that both $`H`$ and $`Z^2`$ return incorrect values for B0531+21 in some VPs where the pulsar is weak.

For each pulsar in Table II, we have searched all VPs for which the pulsar is less than 10o off-axis, with the number searched given in Table II. Using $`f`$ and $`\dot{f}`$ errors as quoted above, we only detect B0531+21 in the last seven of 14 VPs. Widening the range of $`f`$ and $`\dot{f}`$ searched, we detect the pulsar in all VPs. While B0531+21 is strong enough that it is still the most significant signal over the wider range, this would not be the case for the new pulsars. We detect B1706–44 in two of four VPs. If we do not use PSF weighting, only one of these detections is significant. Generally, detection significance for the B0531+21 and B1706–44 trials varies as expected with exposure time and off-axis angle, but we do not detect B1706–44 in one viewing period which has similar exposure and off-axis angle to another in which it is detected with high significance. For comparison, J1706–44 has a a flux just less than that of 3EG J1746–2851 but higher than those of the other unidentified sources listed in Table II.

From two of the Table II pulsars, we detect tantalizing signals. In three of the 14 VPs searched for J1747$``$2958, the highest significance periodicities have similar pulse profile shapes (see Figure 2). This unique shape is not seen in any observations of any other pulsars. While the component amplitude ratio and widths vary for the different VPs, we see similar variations for the B0531+21 detections. It is unclear why this periodicity is not seen in any of the other VPs. The timing model may not be accurate enough to extrapolate to the earliest VPs, but there are some VPs between the first and second detections in which we would expect to see the signal but do not. This may be due to variations in the flux of the source. Similarly, in two of the eight VPs analyzed for J2021+3651, we detect high significance periodicities with similar profile shapes (see Figure 2). Again, why we do not detect the pulsar in other VPs is unclear. While we are working on more thorough analyses of these pulsars, determining whether they are in fact gamma-ray pulsars may have to wait for future observations with a current radio ephemeris. From several other pulsars in Table II we detect marginally significant periodicities. However, none of the profile shapes are repeatable across VPs.

## 5 Future

There are several improvements that we would like to make to our current population model. Obviously, all pulsars will not beam towards us with the same solid angle. While we believe that this assumption is sufficient for making the large-scale predictions of Table III, it does not allow us to explain the detection or non-detection of individual objects. We would also like to incorporate a more realistic pulsar spatial distribution, incorporating spiral arms. This will not make a large difference for the EGRET and AGILE detection statistics but may be important for GLAST, which will probe the population of weak, distant pulsars. Our model currently ignores gamma-ray emission from millisecond pulsars. But with the Kuiper et al. (2000) announcement of a probable detection of pulsed gamma rays from pulsar J0218+4232, this issue must be understood, as it could be important for the detectable populations of AGILE and GLAST. Finally, before the launch of GLAST and AGILE, we must further optimize pulsar search techniques so that we determine the most efficient ways of performing blind searches for gamma-ray pulsars on these data.

## Acknowledgements

We thank Fernando Camilo and Mallory Roberts for providing up-to-date radio ephemerides.


format: text

text:
# Ballistic trajectory: parabola, ellipse, or what?

## I Introduction

The trajectory of a particle moving without drag under the influence of the gravitational field of a perfectly spherical earth is a conic section. Since a cannonball, for example, is on a bound orbit, its flight from muzzle to target is part of an ellipse, as shown on the left in Fig. 1. The center of the earth is at the (distant) focus of the highly eccentric ellipse, and the part of the ellipse relevant to cannonball flight is a small segment near the apogee. Introductory texts, however, treat ballistic motion using Galileo’s approach. Specifically, the horizontal and vertical components of the motion of the projectile are separated, and the trajectory is given by a parabola. The parabola as a limit of an ellipse was noted by Newton in the Principiaprincipia :

> If the ellipsis, by having its centre removed to an infinite distance, degenerates into a parabola, the body will move in this parabola; and the force, now tending to a centre infinitely remote, will become equable. Which is Galileo’s theorem.

Many of the introductory textbooks we have reviewed do not discuss the nature of the approximation, i.e., they do not discuss the conditions under which the flat-earth parabola is an accurate approximation for what is really a central force ellipse. Some textbooks, however, do make explicit statements about the approximation, statements that seem plausible, but turn out to be incorrect. Specifically, we found three classes of conditions which are used in those textbooks that do discuss the validity of the approximation. In the first class it is assumed that the maximal height of the trajectory above the surface of the earth is small compared with the radius of the earth (or, equivalently, that the magnitude of the gravitational acceleration does not change by much along the trajectory)height ; in the second class it is assumed that the range of the trajectory is small compared with the radius of the earth (or, equivalently, that the direction of the gravitational acceleration does not change by much along the trajectory)range ; in the third class it is assumed that both the height and the range are small compared with the radius of the earth (both the magnitude and the direction of the gravitational acceleration do not change by much)height\_and\_range . Other textbooks refer vaguely to neglecting the curvature of the earth in justifying the approximation, but do not quantify the condition vague .

In this paper we show that the correct condition under which the flat-earth picture is a good approximation for the trajectory, is that the maximal curvature of the trajectory is much greater than the curvature of the surface of the earth. This condition coincides with the three classes of naive conditions for many, but not all, trajectories. There are trajectories for which all the three naive classes of conditions are satisfied, but for which the approximation fails to be valid.

To “flatten the earth” without any approximation, we introduce geocentric polar coordinates $`r,\varphi `$ in the plane of the orbit (as shown on the left of Fig. 1) and we plot them (on the right in Fig. 1) as if they were Cartesian coordinates $`x,y`$. We let $`R`$ denote the radius of our perfect earth, and we take $`h`$ to be the height, above the earth, of the peak of the trajectory. The apogee of the particle orbit is then at a distance $`r=R+h`$ from the center of the earth. For convenience, we choose our zero of $`\varphi `$ so that it coincides with the apogee of the orbit.

In order to view the space above the earth in this way, as if the earth were flat, we need to distort the coordinates in somewhat the same manner as a map maker. We choose to do this with the explicit relationship notunique

$$x=(R+h)\varphi y=rR.$$
(1)
and to plot $`x,y`$ as if they were Cartesian coordinates.

In this paper we shall focus on two questions about exact ballistic trajectories. The first is the accuracy of the parabolic approximation. The naive expectation of most students (and many teachers) is that the approximation is valid if the trajectory is “small” compared to the radius of the earth, i.e. , if neither $`h`$ nor the range of the trajectory is comparable to the earth radius $`R`$. We shall see that this naive expectation is too naive.

Our second question involves the vertical nature of gravitation. In the “true” picture of gravitational acceleration, on the left side of Fig. 1, the gravitational force is radially directed towards the center of the earth. These radial lines are converted to vertical lines by the redrawing of the coordinates on the right hand side of Fig. 1, and by the transformation in Eq. (1). It would seem, therefore, that gravitational acceleration must be acting vertically in the picture on the right hand side of Eq. (1). The introductory textbooks add to this the assumption that $`g`$ is constant, and prove that a particle being acted upon only by gravitational forces moves in a parabola. This argument leads unavoidably to the conclusion that the particle trajectory deviates from a parabola only due to the slight variation of $`g`$ with altitude. We shall see that again the naive expectation is not correct.

## II Elliptical orbits in flat-earth coordinates

To describe the orbit of a unit mass particle we start by letting

$$L=r^2\frac{d\varphi }{dt}$$

represent the particle angular momentum per unit particle mass. With standard techniquesanymechbook  we can write the equation of the particle’s orbit, in polar coordinates, as

$$r=\frac{1}{\frac{GM}{L^2}+\left(\frac{1}{R+h}\frac{GM}{L^2}\right)\mathrm{cos}\varphi },$$
(2)
which is the equation of an ellipse (though not in the most familiar form). Here $`M`$ is the mass of the earth and $`G`$ is the universal gravitational constant.

We can convert Eq. (2) to the flat-earth picture simply by introducing the $`x,y`$ coordinates of Eq. (1), and can write the trajectory in the right hand side of Fig. 1 as

$$y+R=\frac{1}{\frac{GM}{L^2}+\left(\frac{1}{R+h}\frac{GM}{L^2}\right)\mathrm{cos}(x/[R+h])}.$$
(3)
To simplify the complex appearance of this equation we introduce two dimensionless parameters. The first is $`ϵh/R`$, a parameter that tells us something about the relative size of the trajectory and the earth. We want our second parameter to be expressed in terms of quantities appropriate to the usual description of a ballistic trajectory. To that end we introduce the parameter $`V_{horiz}`$, the particle velocity at the top of the trajectory, and we combine it with the acceleration of gravity at the surface of the earth $`g`$ to form the dimensionles parameter

$$\alpha \frac{V_{horiz}^2}{gh}.$$
(4)
We next note that $`V_{horiz}=L/(R+h)`$ and $`g=GM/R^2`$, so that $`\alpha `$ is equivalent to $`(L^2/GMh)/(1+ϵ)^2`$. With these equivalences we can rewrite the trajectory in Eq. (3) as

$$\frac{y}{h}=\frac{\alpha (1+ϵ)^2}{1+\left[\alpha ϵ(1+ϵ)1\right]\mathrm{cos}\left(\frac{ϵ(x/h)}{[1+ϵ]}\right)}\frac{1}{ϵ}.$$
(5)
In Eq. (5), the only length scale that explicitly appears is $`h`$, a scale appropriate to the description of ballistic trajectories. But the parameter $`ϵh/R`$ still carries information about the size of the earth. The familiar parabolic trajectory follows when this parameter is taken to be very small. To see this, we can expand the right hand side in powers of $`ϵ`$ and keep only the term to zero order in $`ϵ`$ (keeping $`\alpha `$ fixed). The result is

$$\frac{y}{h}=1\frac{1}{2\alpha }\frac{x^2}{h^2}=1\frac{g}{2V_{horiz}^2}\frac{x^2}{h},$$
(6)
the standard parabolic ballistic trajectory,

We now look more carefully at the exact orbit in the flat-earth coordinates of Eq. (5). In Fig. 2 we present trajectories for $`ϵ=10^4`$, hence for $`h0.63`$ km. The derivation above, of the parabolic approximation, assumes that $`\alpha `$ is of order unity. To see interesting deviations from a parabola, we consider large values of $`\alpha `$ in the figure. For $`\alpha =15,000`$, the deviations from the parabola are striking indeed; the trajectory bends upward. This, of course, is simply the appearance in our coordinates of an elliptical orbit for which the radius of curvature at the apogee is greater than the radius of curvature of the earth’s surface.

What may be most surprising in Fig. 2, is that large deviations from the parabola occur when the range of a trajectory is only around $`90h`$, or about 60 km, a tiny distance compared to the size of the earth. To see why this is so, we can put $`y=0`$ in Eq. (5) and we can solve for $`x`$, half the range of the particle:

$$\frac{x}{h}=\frac{1+ϵ}{ϵ}\mathrm{cos}^1\left(\frac{1ϵ\alpha (1+ϵ)^2}{1ϵ\alpha (1+ϵ)}\right)=\sqrt{2\alpha }\left[1+\frac{ϵ}{2}(\alpha 1)\right]+𝒪\left(ϵ^2\right).$$
(7)
The fractional error in the range of the trajectory, for large $`\alpha `$ is therefore

$$\frac{ϵ}{2}(\alpha 1)\frac{ϵ\alpha }{2}\frac{\mathrm{range}^2}{16Rh},$$
(8)
where we have used the fact that the range is approximately $`2h\sqrt{2\alpha }`$. Equation (8) shows that the deviation from the standard range formula is large when the range is of order $`\sqrt{Rh}`$ or $`R\sqrt{ϵ}`$, and hence can be important for orbits that are much smaller than the scale of the earthloworbit . This conclusion is in good agreement with the $`\alpha =3000,ϵ=10^4`$ example in Fig. 2. For these parameters the parabolic approximation misses the target by $`15\%`$.

There is a nice way of understanding what we have just found. The replacement of the spherical earth surface with a flat surface, in Fig. 1, is justifiable only if the trajectory is more sharply curved than is the earth surface. That is, if the minimal radius of curvature along the trajectory is smaller than the radius of the earth. The radius of curvaturecurvature  $``$ of the trajectory iscurvature\_ellipse

$$:=\frac{\left[1+(dy/dx)^2\right]^{3/2}}{\left|d^2y/dx^2\right|}=\frac{h\alpha }{\left[1+x^2/\left(h\alpha \right)^2\right]^{3/2}},$$
(9)
where we have used Eq. (6). The condition that the maximum curvature of the trajectory (the curvature at $`x=0`$) is much greater than the curvature of the earth’s surface is that $`h\alpha /R1`$. This is violated when $`ϵ\alpha 1`$ for which we have, from Eq. (8), a large deviation from the parabolic range prediction. We can now understand this in terms of radius of curvature of the trajectory. For a small $`h`$, high velocity orbit, a projectile can move “just above” the surface of the earth on a trajectory with a range much larger than its height $`h`$ (but much smaller than $`R`$). The curvature of that trajectory in space is mostly due to the curvature of the earth. In the flat earth picture the true trajectory would seem to have almost constant height, and would greatly deviate from a parabolalatus\_rectum . It should be noted that the meaning of a flat earth trajectory disappears if $`ϵ\alpha >1`$. It follows from Eq. (9) that for $`ϵ\alpha >1`$ the radius of curvature at the apogee exceeds that of the earth. This suggests that the elliptical or hyperbolic trajectory of the projectile will not intersect the earth surface. That this is, in fact, true can be verified by looking for solutions of $`y=0`$ in Eq. (5); for $`ϵ\alpha >1`$, there are no solutions with positive $`ϵ`$.

There is also an interesting non-geometrical way of viewing the condition $`\alpha ϵ1`$. From the definitions of $`\alpha `$, $`ϵ`$, and $`g`$, we have

$$\alpha ϵ=\frac{V_{horiz}^2}{gR}=\frac{V_{horiz}^2}{V_{circ}^2}=\frac{2V_{horiz}^2}{V_{esc}^2},$$
(10)
where $`V_{circ}`$ is the velocity of a circular orbit with radius $`R`$ (i.e. , just above the earth surface), and $`V_{esc}`$ is the escape velocity from the earth surface. The condition $`\alpha ϵ1`$, then, is simply the condition that the motion is very slow compared to typical “orbital” (as opposed to “trajectory”) motions. For trajectories with $`V_{horiz}`$ comparable to $`V_{esc}`$ we should not be surprised that the deviation from a parabolic projectile trajectory is significant. We may, however, be surprised that the range of that trajectory may be short, much less than the earth radius. The key idea here is that the range is kept short by imposing a small value of $`h`$. For sufficiently small values of $`h`$ the high velocity projectiles hit the earth’s surface before they get a chance to show how far they would have gone in their large elliptical orbits. In order for this to happen, in order for the range to be kept small, the value of $`h`$ must be reduced as $`V_{horizon}`$ is increased. The limiting case is $`V_{horizon}=V_{esc}/\sqrt{2}`$, the maximum for which the trajectory has a range (i.e., for which the trajectory intercepts the earth’s surface). As this limit is approached, the condition on $`h`$, for the range to be much less than the earth radius, is $`ϵ\left[12(V_{horizon}/V_{esc})^2\right]`$.

## III Is gravity vertical in the flat-earth picture?

In Section I we argued that gravitational acceleration must be vertical in the flat earth picture. It is apparent in Fig. 2 that this cannot be correct. For the $`\alpha =15,000`$ trajectory, the particle is accelerating upward. And this particle is at a height of only a kilometer or so above the earth’s surface, where the acceleration of gravity is certainly close to 9.8 m/sec2 and is even more certainly downward.

Another clear indication that “vertical acceleration” is not the whole story, is the fact that the horizontal velocity $`dx/dt`$ changes in time. Figure 3 shows the small, but nonzero, fractional change in the horizontal velocity along the orbit.

To understand the results in this figure we note that

$$\frac{1}{V_{horiz}}\frac{dx}{dt}=\frac{(R+h)}{V_{horiz}}\frac{d\varphi }{dt}=\frac{(R+h)L}{V_{horiz}r^2}=\frac{(R+h)^2}{r^2}.$$
(11)
In this simple calculation it is easy to see why the horizontal velocity changes: it is proportional to $`L/r^2`$. The angular momentum $`L`$ is constant, but $`r`$ is not constant except in the case of the circular orbit for which $`\alpha =\alpha _{crit}R^2/h(R+h)`$. As the particle moves to smaller radius (in the case $`\alpha <\alpha _{crit}`$) the horizontal velocity must increase; as the particle moves to larger radius (in the case $`\alpha >\alpha _{crit}`$) the particle velocity must decrease. In this mathematics it is easy to see why the horizontal velocity changes during the motion.

At a deeper level, the failure of naive intuition is due to the tacit expectation that the equations of motion for a unit mass particle can be put in a form

$$\frac{d^2x}{dt^2}=F_x\frac{d^2y}{dt^2}=F_y,$$
(12)
where $`F_x`$ and $`F_y`$ are only functions of position $`x,y`$. Since particles released from rest will fall vertically downward we conclude that $`F_x=0`$, and gravity is vertical.

We are misled to this conclusion by the Cartesian appearance of the $`x,y`$ coordinates on the right side of Fig. 1. The actual equations of motioneqsofmotion , turn out to be

$`{\displaystyle \frac{d^2x}{dt^2}}`$ $`=`$ $`{\displaystyle \frac{2}{(y+R)}}{\displaystyle \frac{dx}{dt}}{\displaystyle \frac{dy}{dt}}`$ (13)
$`{\displaystyle \frac{d^2y}{dt^2}}`$ $`=`$ $`{\displaystyle \frac{GM}{(y+R)^2}}+{\displaystyle \frac{(y+R)}{(R+h)^2}}\left({\displaystyle \frac{dx}{dt}}\right)^2.`$ (14)
In the case of a vertical trajectory ($`dx/dt=0`$) the equations are in complete agreement with our expectations. When $`dx/dt0`$, however, the velocity dependent terms, absent in Eqs. (12), are responsible for the failure of intuitonnotG ; neglect\_vel . These velocity dependent terms are due to the fact that the $`x,y`$ coordinate lines are actually curved lines through physical space. We can draw them as straight lines, but the way in which we draw them, or the symbols we use to represent them, do not change the fact that they are curved in physical space. This curvature of spatial coordinates comes with its usual kinematical consequences. The $`(dx/dt)^2`$ in Eq. (14), for example, is just the usual centripetal term disguised by unusual coordinate names.

## IV Conclusions

It is plausible to expect that the uniform-gravity, flat-earth gravity parabola is an accurate approximation if the height $`h`$ of the trajectory, and the range of the trajectory, are both much less than the earth radius $`R`$. In terms of our parameterization these conditions are, respectively $`ϵ1`$ and $`\alpha ϵ^21`$. If we are to have the kinematics correct, as well as the shape of the orbit, then clearly the first condition is necessary; if $`h`$ is not small compared to $`R`$ gravity cannot be of uniform strength.

The second condition, however, is not sufficiently strict. The correct second condition for the validity of the flat-earth approximation is that the maximum curvature of the trajectory be much greater than the curvature of the earth’s surface. In our notation this corresponds to orbits that satisfy $`\alpha ϵ1`$. This condition is not satisfied for a class of trajectories that do satisfy $`ϵ1`$ and $`\alpha ϵ^21`$. That class consists of low height, high velocity motions. The importance of the high velocity feature of these orbits can be seen in the fact that our condition $`\alpha ϵ1`$ is equivalent to the condition that the velocity of the orbit is much less the escape velocity from the earth surface.

## V Acknowledgment

We thank Benjamin Bromley for pointing out the low-velocity interpretation, given in Eq. (10), for the parabola condition. This work has been partially supported by the National Science Foundation under grant PHY0244605.


text:
# Slow dynamics of a confined supercooled binary mixture II: 𝑄 space analysis

## I Introduction

Modifications of vitrification properties of liquids under confinement are of interest for relevant technological and biophysical problems pablo ; grenoble . While experiments show that the glass transition is not usually suppressed when the glass former is confined, changes in the temperature of the transition and in dynamical properties of liquids upon supercooling are induced by both the geometrical confinement and the interaction with the substrate grenoble .

In the study of the glass transition phenomenology the mode coupling theory (MCT) of the evolution of glassy dynamics plays a fundamental role. MCT is able to describe the slowing down of supercooled liquids in a region of mild supercooling on approaching a temperature $`T_C`$ gotze1 ; gotze2 ; gotze3 . This temperature marks a crossover from a region where dynamics is mastered by the transient caging of nearest neighbor to a lower temperature region where hopping processes dominate. In the ideal version of the theory, where hopping is not included, the non linear feedback mechanisms in the microscopic dynamics of the particles become so strong that they lead to the structural arrest of the system. Close to $`T_C`$ MCT asymptotic solutions show analytic behaviors of the density correlator from which parameters of the theory can be extracted. These behaviors have been effectively found in experiments and simulations on supercooled liquids in regions close to $`T_C`$ where hopping was negligible gotze3 .

The role of the cooperative dynamics of the particles in the glass transition is still under discussion and there is an increasing activity on the definition and observation of dynamical heterogeneities adam-gibbs . Cooperativity in the dynamical behavior of the particles increases as the system is approaching the glass transition, suggesting the idea of the existence of a dynamic correlation length donati1 ; donati2 ; doliwa1 ; doliwa2 ; berthier .

Molecular dynamics studies of model liquids in restricted geometries intended to assess the applicability of MCT can give an important contribution to the characterization of vitrification processes in confinement gallo1 ; gallo2 ; baschnagel ; kobconf1 ; kobconf2 .

Liquids confined in network of interconnected pores with a large value of porosity, as Silica Xerogels, can be appropriately studied with models where the confining solid is built as a disordered array of frozen microspheres rosinberg ; monson ; page . We consider here one of such system. A liquid Lennard Jones binary mixture (LJBM) composed of $`80\%`$ of particles A and $`20\%`$ of smaller particles B embedded in an off lattice matrix of soft spheres. We performed Molecular Dynamics (MD) simulations of the confined LJBM upon supercooling in order to test the predictions of MCT.

The bulk phase of this LJBM behaves à la Mode Coupling kob ; kob1 ; kob2 . A preliminary analysis on our system euro  showed evidences that predictions of MCT hold also for the confined LJBM. A detailed study on the dynamics carried out in the direct space has shown two important differences with respect to the bulk upon supercooling noi . The smaller B particles tend to avoid the soft sphere interfaces on lowering temperature and correspondingly their diffusion coefficient becomes lower than that of the A particles, at variance with the bulk where there is no such an inversion. Hopping processes are markedly present for both A and B particles already above $`T_C`$ while in the bulk system only B particles show consistent hopping.

As a possible consequence of these differences, we observed that the range of validity of the power law fit to the diffusion coefficient as predicted by MCT in the long time region, late $`\alpha `$ relaxation region, appears much more limited than in the bulk.

A more refined assessment of the MCT validity can be obtained by a quantitative comparison between the asymptotic predictions of the theory and the results of our computer simulation in the $`(Q,t)`$ space, where the density correlators of a liquid testing MCT display clear signatures of the approach to $`T_C`$. In this space important scaling relations of the theory can be verified. For these reasons we carry on here an analysis of the behavior of the self intermediate scattering functions, $`F_S(Q,t)`$.

In the next section we recall the predictions of the idealized version of MCT connected to our analysis. In Sec. III we describe our model system and Molecular Dynamics details. In Sec. IV we perform a test of MCT properties on our system. In particular we analyze the static structure factor and the self intermediate scattering function. For this latter quantity we perform numerical fit to MCT analytical laws and extract parameters of the theory. Last section is devoted to conclusions.

## II MCT scenario for the evolution of the structural relaxation

The MCT for the evolution of glassy dynamics developed originally to deal with the cage effect provides a model for an ideal liquid-glass transition. The evolution of the structural relaxation is described in terms of the asymptotic behavior of density correlators near singularities which take place close to a cross-over temperature $`T_C`$.

For a density correlator or any other correlator that couples with density fluctuations, $`\varphi (Q,t)`$, MCT predicts upon supercooling a diversification of relaxation times related to the nearest neighbor caging. In the short time regime the correlator show a fast relaxation associated with the ballistic motion of the particles. The high temperature correlator relaxes then exponentially for longer times. On lowering temperature for intermediate times a shoulder appears in the relaxation law. As the temperature is further lowered the shoulder enhances and becomes a plateau and the correlation function clearly shows a two step relaxation. The presence of the plateau is related to the trapping of the particle in the cage formed by its nearest neighbors. The time interval in which the correlation function is close-to and above the plateau is called the $`\beta `$ relaxation region. The region that starts when the correlation function leaves the plateau is called the $`\alpha `$ relaxation region.

The long time limit of $`\varphi (Q,t)`$ is the non ergodicity parameter $`f_Q(T)`$. This parameter exhibits a singularity at a temperature $`T_C`$ where the system undergoes an ideal transition from an ergodic to a non ergodic behavior. $`f_Q(T)`$ jumps discontinuously from $`0`$ above $`T_C`$ to $`f_Q^c`$ at $`T_C`$. On approaching $`T_C`$ from below

$$f_Qf_Q^c\sqrt{|ϵ|}$$
(1)
where $`ϵ=(TT_C)/T_C`$ is the small parameter of the theory and $`f_Q^c`$ is called the critical non ergodicity parameter. On approaching $`T_C`$ from above an effective $`f_Q`$ can be identified with the height of the plateau of the correlator. In this range of temperatures the values of $`f_Q`$ are very close to $`f_Q^c`$ and show a very mild temperature dependence

$$f_Qf_Q^cO(ϵ)$$
(2)
In the asymptotic limit where $`ϵ<<1`$ MCT predicts a factorization

$$\varphi (Q,t)f_Q^c=h_QG(t)$$
(3)
where $`h_Q`$ is a positive factor which does not depend on time. The so called $`\beta `$ correlator, G(t), obeys the scaling relation

$$G(t)\sqrt{|ϵ|}g(t/t_ϵ)$$
(4)
where the time scale $`t_ϵ`$ has been introduced

$$t_ϵ\frac{t_0}{|ϵ|^{1/2a}}.$$
(5)
Here $`t_0`$ is a microscopic characteristic time of the system and the exponent $`a`$ is called the critical exponent. In the time scale which corresponds to the decay of the correlator towards the plateau, $`tt_ϵ`$

$$g\left(t/t_ϵ\right)=\left(t_ϵ/t\right)^a$$
(6)
For the correlator decay of the liquid from the plateau

$$g\left(t/t_ϵ\right)\left(t/t_ϵ\right)^b$$
(7)
which is valid for times larger than the time scale of the $`\beta `$ relaxation and much smaller than the $`\alpha `$-relaxation time $`\tau `$, $`t_ϵt<<\tau `$. The exponent $`b`$ is known as the von Schweidler exponent. The critical and the von Schweidler exponents are related by the formula

$$\lambda =\frac{[\mathrm{\Gamma }\left(1a\right)]^2}{\mathrm{\Gamma }\left(12a\right)}=\frac{[\mathrm{\Gamma }\left(1+b\right)]^2}{\mathrm{\Gamma }\left(1+2b\right)}$$
(8)
where $`\lambda `$ is called the exponent parameter and $`\mathrm{\Gamma }(x)`$ is the Euler Gamma function.

In the late part of the $`\alpha `$-relaxation regime the shape of the correlator of a liquid approaching the MCT crossover temperature gotze2  can be represented by means of the Kohlrausch-Williams-Watts (KWW) analytical function

$$\varphi (Q,t)=f_Q\mathrm{exp}((t/\tau )^\beta )$$
(9)
where $`\beta `$ is the Kohlrausch exponent.

The $`\alpha `$-relaxation time $`\tau `$ is expected from MCT to follow the power law

$$\tau ^1(TT_C)^\gamma $$
(10)
Both the exponent $`\gamma `$ and the crossover temperature $`T_C`$ should coincide with the ones appearing in the temperature dependence of the diffusion coefficient $`D(TT_C)^\gamma `$.

An important prediction of MCT is the validity of the time-temperature superposition principle. It states that in the low temperature regime, approaching from above the MCT crossover temperature, it is possible to time-rescale the correlators evaluated at different T into a single master curve:

$$\varphi (t)=\stackrel{~}{\varphi }(t/\tau (T))$$
(11)
$`\stackrel{~}{\varphi }`$ is the master function which in the late $`\alpha `$-relaxation limit becomes the KWW function (9).

In the late $`\beta `$, early $`\alpha `$ relaxation regime Eq. 3 can be written as

$$\varphi (Q,t)=f_Q^ch_Q(t/\tau )^b$$
(12)
This equation is called the von Schweidler law and $`\tau `$ indicates a characteristic time of the system.

The exponent $`\gamma `$ in Eq. 10 is related to the exponents $`a`$ and $`b`$ by the equation

$$\gamma =\frac{1}{2a}+\frac{1}{2b}$$
(13)
While both $`\tau `$ and the exponent $`\beta `$ are $`Q`$ dependent, the exponent $`b`$ should not depend on $`Q`$. These two exponents are related: $`\beta >b`$ and $`\beta (q\mathrm{})=b`$ fucs .

## III Model and Molecular Dynamics

The model studied consists in a rigid off lattice matrix of 16 soft spheres where a liquid $`A_{80}B_{20}`$ Lennard Jones binary mixture (LJBM) kob1 ; kob2  is embedded. The mixture is composed of 1000 particles. In the following LJ units will be used. The parameters of the LJBM are $`ϵ_{AA}=1`$, $`\sigma _{AA}=1`$, $`ϵ_{BB}=0.5`$, $`\sigma _{BB}=0.88`$, $`ϵ_{AB}=1.5`$, and $`\sigma _{AB}=0.8`$. The parameters of the soft sphere potential are $`ϵ_{SA}=0.32`$, $`\sigma _{SA}=3`$, $`ϵ_{SB}=0.22`$, $`\sigma _{SB}=2.94`$. The box length is L=12.6. We performed Molecular dynamics simulations in the NVE ensemble. We studied a set of temperatures ranging from T=5.0 down to T=0.37. The timestep for the lowest temperature is 0.02 and the largest production run consists in $`14`$x$`10^6`$ timesteps. For all temperatures the time range investigated by MD was sufficient for the correlation functions to decay to zero in the long time limit. Further MD simulation details are reported in ref.s noi ; euro . We did not find any qualitative differences in the behavior of the thermodynamic quantities with respect to the bulk upon supercooling noi .

## IV Analysis of the dynamical behavior in $`(Q,t)`$ space

### IV.1 Static structure factor

In the framework of MCT the static structure factor is the input information for MCT equations gotze1 ; gotze2 . This function is expected to depend smoothly on wavevectors and control parameters and it is also expected to show no singularities upon supercooling. In Fig. 1 we report the structure factors $`S(Q)`$ of the confined mixture for the temperatures investigated. They behave consistently with MCT as the functions smoothly vary upon supercooling implying no divergence of correlation length. No signature of phase separation of the mixture are evident in confinement as we already deduced from the analysis of the $`g_{AB}(r)`$ noi . We observe that there are no relevant discrepancies in the confined $`S_{ij}(Q)`$ with respect to the bulk kob2  neither in the shape nor in the peak positions. For $`Q`$ values around the first structure factor peak MCT features are best evident in the density correlators. The peak positions in our confined LJBM show a very mild temperature dependence, especially for low temperatures. We will therefore consider them constant. The values $`Q_{MAX}^A=7.06`$ $`Q_{MAX}^B=5.90`$ will be used to test MCT as a function of temperature in the following. In the case of the bulk the values used were respectively $`Q_{MAX}^{AA}=7.25`$ and $`Q_{MAX}^{BB}=5.75`$.

### IV.2 Self intermediate scattering function

We consider now the self intermediate scattering function (SISF) defined as

$$F_S(Q,t)=\frac{1}{N}<\underset{i}{}e^{i\stackrel{}{Q}\stackrel{}{R}_i(0)}e^{i\stackrel{}{Q}\stackrel{}{R}_i(t)}>$$
(14)
This correlator can be used to test the predictions of MCT close to $`T_C`$ as described in Sec. II. In Fig. 2 the SISF of our confined system is plotted for all considered temperatures and for both particle types at the peak of the structure factors. At high temperatures the correlators relax exponentially for longer times. On lowering temperature we observe the onset of the plateau related to the two step relaxation typical of a supercooled liquid.

An onset temperature for the plateau of $`T0.6`$ can be identified for our system. The corresponding temperature for the bulk is $`T1.0`$. Analogous to bulk supercooled liquids for the lowest temperature we observe that correlators stretch several decades in time. We also see that A and B particles relax in a similar fashion but larger A particles relax faster at this wavevectors, especially for lower temperatures. This counterintuitive phenomenon that happens only in confinement is possibly connected to packing constraints. We have in fact detected that on lowering temperature smaller B particles tend to avoid the soft sphere interfaces and preferentially move in the center of the interstices. Correspondingly diffusion coefficients of A and B particles extracted from the mean square displacement (MSD) show an inversion as supercooling progresses noi .

### IV.3 Kohlrausch-Williams-Watts analysis of the self intermediate scattering function

In Fig. 2 together with the SISF we show the fit to the KWW function, as defined in Eq. 9, performed for temperatures $`T<0.8`$ in the $`\alpha `$ relaxation region for both A and B particles. A remarkable agreement with the data is found. From the fit the KWW relaxation time $`\tau `$, the stretching exponent $`\beta `$ and the effective non ergodicity parameter $`f_Q`$ are extracted.

MCT predicts both $`D`$ and $`\tau ^1`$ to follow the same power law, see Eq.10. The relaxation times $`\tau `$ are shown in the main frame of Fig. 3. In the same picture we show the diffusion coefficients $`D`$ previously reported and obtained from the slope of the MSD noi . We also display in the figure the fit to the MCT predicted behavior. Values of $`T_C`$ and $`\gamma `$ obtained from the fit are reported in Tab.1. We observe that both in bulk and confined system the crossover temperature $`T_C`$ extracted from $`D`$ and $`\tau `$ and from A and B particles coincide, $`T_C0.36`$ A substantial reduction of $`T_C`$ is observed in confinement. The bulk value reported in literature is $`T_C0.43`$.

The exponents $`\gamma `$ which determine the slope of the fit are the same for A and B particles both for $`\tau `$ and $`D`$. This agrees with the prediction of MCT that these exponents should be universal and do not depend on the species of the particles involved. A discrepancy with respect to MCT prediction is represented by the different value of $`\gamma `$ between $`D`$ and $`\tau `$. As we can see from the table the same discrepancy affects the mixture also in the bulk phase and it may be due to long time activated processes causing a breakdown of the predictions of the theory only for quantities calculated at extremely long times noi . In fact in confinement, where hopping is consistent for both A and B particles, a reduction of circa $`35\%`$ is observed for $`\gamma `$ in going from $`\tau `$ to $`D`$. In the bulk phase, where hopping is more consistent for B particles, the reduction observed is $`30\%`$ for B and $`20\%`$ for A particles.

A substantial shrinkage of the range of validity of the theory is found in confinement. In Fig. 3 we in fact see that the MCT predicted power law behavior is verified for $`\tau `$ in the temperature window $`0.58T0.41`$, the same found for the $`D`$ coefficients noi , which corresponds to $`0.153<ϵ<0.631`$ while in the bulk $`0.07<ϵ<1.30`$ kob .

The values of the parameters $`f_Q`$ and $`\beta `$ extracted from the fit are plotted as a function of temperature in the inset of Fig. 3. We note that the effective non ergodicity parameter $`f_Q`$ only slightly depends on temperature in agreement with MCT that predicts a consistent temperature dependence of this factor only below $`T_C`$. The stretching parameter $`\beta `$ decreases upon lowering the temperature as also expected, but we observe no flattening of this exponent on approaching $`T_C`$ as found for bulk systems. Besides, correlators appear much more stretched in confinement. In the bulk $`\beta >0.78`$ while here the lowest values are slightly below $`0.4`$.

We report in Fig. 4 the temporal evolution of the correlators for different wavevectors $`Q`$ ranging from $`2.5`$ to $`14.5`$ for the temperature $`T=0.41`$. In the same picture we also show the fit to the KWW function. A remarkable agreement is found. The height of the plateau strongly depends on $`Q`$. The same analysis was performed for $`T=0.58,0.48,0.43,0.39`$ (not shown).

The values of the parameters of the fit are reported for all the above temperatures in Fig. 5. All the parameters show a monotonic behavior in $`Q`$. In particular at high temperature the relaxation time $`\tau `$ spans two decades. It is worth noting that the $`\beta `$ exponent is converging, for $`TT_C`$ and for $`Q\mathrm{}`$ to the value of $`b`$ for both types of particles as expected fucs . Finally we note that the prefactor $`f_Q`$ is slightly $`T`$ dependent for all $`Q`$.

In Fig. 6 the $`\tau `$ values extracted from the fit are reported as a function of the shifted temperature $`TT_C`$. In the temperature range where the power law behavior is effective, that for our system is $`0.41T0.58`$, $`\gamma `$ is expected to be independent on $`Q`$. In Fig. 6 we effectively observe that the curves are parallel to each other and parallel to the straight line overlayed that corresponds to the value of $`\gamma =2.9`$ extracted from the fit performed at the peak of the structure factor (see Tab.1).

### IV.4 von Schweidler test for the self intermediate scattering function

According to Eq. 11 the time dependence of the correlation functions can be rescaled by the relaxation time $`\tau `$ extracted from the fit to the KWW. The curves for the lowest temperatures are expected to collapse onto a mastercurve in the late $`\beta `$, early $`\alpha `$ relaxation region. In Fig. 7 we report for both particle types the master plot of the SISF of Fig. 2. In the same figure we show the fit to the von Schweidler functional form of Eq. 12.

The values of the parameters obtained for A and B particles from the fit are:

$$f_{Q,A}^c=0.72,h_{Q,A}=0.35,b_A=0.355$$

$$f_{Q,B}^c=0.79,h_{Q,B}=0.44,b_B=0.350$$

In the bulk the values of the exponent $`b`$ are respectively $`0.51`$ for A and $`0.46`$ for B particles while we observe that in confinement $`b`$ values are lower and almost the same. This might indicate that hopping is influencing also the earlier $`\alpha `$ relaxation region.

An important prediction of MCT is that the von Schweidler $`b`$ value should be independent of the wave vector $`Q`$. This test was performed in the bulk kob2 , where the authors stressed a small dependence of $`b`$ between $`Q_{MAX}`$ and $`Q_{MIN}`$, where $`Q_{MIN}`$ is the position of the first minimum after the peak of the static structure factor. In Fig. 8 the values for $`f_Q^c`$ and $`b`$ are presented as a function of $`Q`$ for A particles. The values for $`f_Q`$ and $`\beta `$ for the KWW and the $`S(Q)`$ are also shown. The $`f_Q`$ factors are in perfect agreement with the Gaussian behavior predicted by the MCT and found in the bulk simulations. We can therefore consider the $`f_Q`$ KWW factors as a good approximation of the critical non-ergodicity parameter, $`f_Q^c`$, at such low temperature. The $`b`$ parameter decreases monotonically, except between the first $`S(Q)`$ maximum and the following minimum where it is practically constant, similar to the bulk. We note that the $`\beta `$ values are always higher than the relative $`b`$ values, this is a consequence of the fact that the KWW can describe the late $`\alpha `$ but not the late $`\beta `$, early $`\alpha `$ regime where the data are well fitted by the von Schweidler law.

In our correlators oscillations, possibly due to finite size effects, in the time zone where the power law of Eq.6 holds, prevents us from calculating the exponent $`a`$ of MCT directly. It is still possible to evaluate the parameter $`a`$ from the closure relation of Eq. 13. Using the $`\gamma `$ calculated from $`\tau `$ and the von Schweidler exponent $`b`$ found previously, the values of the critical exponent $`a`$ for the two species are respectively

$$a_A=0.335,a_B=0.342$$

These values are quite high, although within the topmost value indicated by the theory $`a0.395`$.

## V Summary and Conclusions

We have studied the dynamics of a model glass former in confinement. This study is intended to assess the validity of MCT for supercooled liquids embedded in strongly repulsive and highly confining environments. Our model consists in a LJBM embedded in a disordered matrix of soft spheres. It represents a possible modellization for liquids hosted in Silica Xerogels.

The analysis performed in $`r`$ space and previously reported noi  showed that the MSD behavior can be interpreted for both A and B particles in terms of the onset of the cage effect and that the diffusion coefficients $`D`$ extracted from the MSD can be fitted with the power law predicted by MCT. There is however a clear evidence that the range of validity of the MCT predictions in the late $`\alpha `$ relaxation region suffers a reduction of $`60\%`$ with respect to the bulk. In terms of the small parameter of the theory $`ϵ`$ the bulk range is $`0.07<ϵ<1.30`$ while in the confined system it results to be $`0.153<ϵ<0.631`$. We inferred that the upper bound reduction might be due to the fact that B particles tend to avoid the soft sphere interfaces in lowering temperature. The MCT power law starts in fact to hold for $`D`$ only when most B particles are avoiding the soft spheres noi . The lower bound reduction is to be connected to the presence of important hopping processes more marked with respect to the bulk, which tend to introduce deviations from the ideal MCT predictions earlier than in bulk.

We have presented here an analysis in $`(Q,t)`$ space of the confined LJBM by considering the SISF and we have tested the main scaling relations of MCT.

The static structure factor is similar to that of the bulk phase. The SISF is well fitted with the stretched exponential function and the power law behavior of the $`\alpha `$-relaxation times $`\tau `$ extracted for both A and B particles is consistent with MCT. Power law fit to both $`D`$ and $`\tau `$ as a function of temperature give the same value for the crossover temperature $`T_C=0.355`$. We observe a reduction of circa $`20\%`$ of $`T_C`$ in going from bulk to confined which can be ascribed to the existence of an upper bound in confinement for the size of the domains of cooperative dynamics. A discrepancy observed in the values of the exponents $`\gamma `$ obtained for $`D`$ and $`\tau `$ is to be connected to the more important presence of hopping processes in confinement. Hopping processes influence in fact the behavior of $`D`$ more than that of $`\tau `$ causing a decrease of the value of the exponent. This discrepancy was present also in the bulk.

The existence of a confining disordered host structure seems to be connected also to the lowering of the stretching parameter possibly caused by a larger distribution of relaxation times with respect to the bulk.

The behavior of the SISF as a function of $`Q`$ is very similar to that of the bulk. All parameters show a monotonic behavior as a function of $`Q`$ and the stretching exponent approaches $`b`$ for high $`Q`$. The exponent $`\gamma `$ shows no $`Q`$ dependence.

The SISF plotted against $`t/\tau `$ shows a scaling behavior. The mastercurve can be fitted with the von Schweidler function in the late $`\beta `$, early $`\alpha `$ region. Analogous to the bulk, the exponent $`b`$ is practically constant between the first maximum and the following minimum of the structure factor. Similar to $`\beta `$ a substantial reduction of $`b`$ is observed in going from bulk to confined LJBM.

The exponent parameter $`\lambda `$ of the system can be calculated from Eq. 8. Contrary to MCT prediction we get different values of $`\lambda `$ on using the exponent $`a`$ or the exponent $`b`$ in Eq. 8

$$\lambda _{a,A}=0.680,\lambda _{b,A}=0.871$$

$$\lambda _{a,B}=0.663,\lambda _{b,B}=0.873$$

This inconsistency is present also in the bulk but it is much less severe. The value of $`\lambda `$ reported for the bulk is $`0.79\pm 0.04`$. Here we have a much larger uncertainty that leads to $`0.8\pm 0.1`$.

Considering the strong confinement experienced by the LJBM, the agreement between MCT predictions and the behavior of our confined system as a function of $`T`$ and $`Q`$ is most remarkable. Differences in behavior with respect to the bulk on approaching $`T_C`$ are in fact to be ascribed to the presence of more relevant hopping effects.

The possibility to connect the decrease of $`T_C`$ and the enhancement of hopping effects in confinement with a key role played by dynamical heterogeneities in the glass transition scenario is an issue that would deserve to be addressed in future studies.


text:
# Quasiparticle Relaxation Dynamics in Two Distinct Gap Structures: La0.67Ca0.33MnO3 and LaMnO3

## Abstract

The spin and quasiparticle relaxation dynamics in La0.67Ca0.33MnO3 and LaMnO3 single crystals and thin films is investigated as a function of temperature and magnetic field by pump-probe spectroscopy. The unusually slow spin-lattice relaxation dynamics is governed by the temperature- and magnetic field- dependent pseudogap in La0.67Ca0.33MnO3, and by the temperature-independent Jahn-Teller gap in LaMnO3. Our results show that the coupled dynamics of charge, spin and lattice is strongly correlated with the distinct gap structures in these manganites.

Perovskite manganites, R1-xAxMnO3 (R and A are rare- and alkaline-earth ions, respectively) have attracted much attention because of the large negative magnetoresistance and the magnetic-field-induced phenomena. The strong correlation between spin, charge, orbital, and lattice degrees of freedom leads to complex phase diagrams and to the coexistence of various forms of ordering. Of particular interest in this context has been the observation of a temperature-dependent pseudogap in a variety of colossal magnetoresistance (CMR) oxides and its relationship to the ferromagnetic (FM) metal-to-insulator (MI) transition . Key to the pseudogap and many other properties of the manganites appear to be the existence of nanoscale charge/orbital fluctuations, which cooperate with Jahn-Teller distortions and compete with the electron itineracy favored by double exchange . A profound knowledge of these ordering processes and their mutual interactions is essential for a better understanding of CMR. Dynamical information on the electronic states of the pseudogap can be obtained, if charge or spin excitations are introduced via pulsed photoexcitation, and the relaxation dynamics of the quasiparticles (QP’s) is studied on a picosecond time scale. Ultrafast optical pump-probe spectroscopy has provided significant insight in the dynamics of spin, electron and lattice in metals , and recently in transition metal oxides .

Here, we report the time-resolved pump-probe reflectivity measurements of photoexcited spin and QP relaxation dynamics in La0.67Ca0.33MnO3 (LCMO) and LaMnO3 (LMO) single crystals and thin films as a function of temperature and magnetic field. Our results are fourfold: First, both manganites reveal an unusually slow ($``$ 10 $`\mu `$s) carrier relaxation process which disappears as the transition temperature is approached from below. This slow decay process is attributed to spin-lattice relaxation of carriers in localized (intragap) states. Second, the quasiparticle relaxation rate is governed by the temperature-dependent pseudogap in LCMO and by the temperature-independent Jahn-Teller distortion gap in LMO. Third, the magnetic-field-induced increase of spectral weight in the pseudogap is evident in both the quasiparticle and spin relaxation dynamics in LCMO. Fourth, the polaron relaxation time, $`\tau _B`$, is temperature-independent in LMO, whereas $`\tau _B`$ follows a $`(T_CT)^\beta `$ dependence in LCMO, indicating that the quasiparticle relaxation dynamics involves fully spin-aligned pseudogap states near the Fermi energy $`E_F`$. Our results show for the first time that the coupled dynamics of charge, spin and lattice is strongly correlated with the distinct gap structures in these manganites.

LCMO and LMO single crystals and thin films were grown by the floating zone method and pulsed laser deposition, the details of which are reported elsewhere . The samples were characterized by electrical resistivity and magnetization measurements. The LCMO single crystal and 400-nm thin film have a Curie temperature $`T_C`$ = 225 K and 260 K, respectively. The LMO single crystal shows a Neel temperature $`T_N`$ = 145 K. For the transient reflectivity measurements the samples were mounted in an optical cryostat. The laser system consists of a Ti:sapphire regenerative amplifier (Spitfire, Spectra-Physics) and an optical parametric amplifier (OPA-800C, Spectra-Physics) delivering 100-fs short pulses at a 1-kHz repetition rate tunable from 600 nm to 10 $`\mu `$m. A two-color pump-probe setup is employed with the pump beam power $`<`$ 6 mW and the probe beam power $`<`$ 1 mW. The unfocused pump beam, spot-diameter $``$3 mm2, and the time-delayed probe beam are overlapped on the sample with their polarization perpendicular to each other. The reflected probe beam is detected with a photodiode detector. A SR250 gated integrator & boxcar averager, and a lock-in amplifier are used to measure the transient reflectivity change $`\mathrm{\Delta }R`$ of the probe beam. For the magnetic-field dependent measurement we used a 9-Tesla superconducting magnetic cryostat from Oxford Instruments, Inc.

Figures 1 (a) and 1 (b) show the time evolution of $`\mathrm{\Delta }R`$ from LCMO and LMO single crystals at different temperatures. The pump and probe wavelength is 800 nm. At low temperature (T $`<`$$`<`$ TC, TN), $`\mathrm{\Delta }R`$ shows initially a fast biexponential decay with relaxation times $`\tau _A`$= 0.5 $``$ 4 ps and $`\tau _B`$ = 50 $``$ 100 ps. The fast process reveals the thermalization of photoexcited quasiparticles which occurs on a time scale $`\tau _{QP}`$ = 0.3 $``$ 3 ps . The second process is characteristic for polaron relaxation .

In addition to the initial fast relaxation processes, a very long-lived negative $`\mathrm{\Delta }R`$ signal remains sufficiently long, decay time $`\tau _{SL}`$ 10 $`\mu `$s, that a negative $`\mathrm{\Delta }R`$ is clearly observable even after 1 ms (Figs. 1(a) and 1(b)). This component is present in both the metallic and insulating phase as observed in the transient optical reflectivity and transmission measurements from charge-ordered Pr0.67Ca0.33MnO3 (PCMO), La0.67Sr0.33MnO3 (LSMO), and Nd0.67Sr0.33MnO3 (NSMO). However, there is no evidence of this component in the transient reflectivity change of the paramagnetic phase of LSMO, which shows a phase transition from ferromagnetic metal to paramagnetic metal at 325 K. Therefore, the long-lived negative component of $`\mathrm{\Delta }R`$ is ascribed to a slow spin relaxation process of the magnetically ordered phase.

Figure 2 illustrates the temperature dependence of $`\mathrm{\Delta }R`$ measured at a time-delay $`\mathrm{\Delta }t=500`$ ps, referred to as $`\mathrm{\Delta }R^{}`$ in LCMO and LMO. The temperature dependence of $`\mathrm{\Delta }R^{}`$, for LCMO single crystal ($`\lambda _{pump}=800`$ $`nm/\lambda _{probe}=800`$ $`nm`$) and thin film ($`\lambda _{pump}=800`$ $`nm/\lambda _{probe}=`$ 5.2 $`\mu m`$) are qualitatively the same and independent of wavelength (Fig. 2a). $`\mathrm{\Delta }R^{}`$ increases with increasing temperature followed by an abrupt drop to zero around $`T_C`$. In contrast, for LMO ($`\lambda _{pump}=800`$ $`nm/\lambda _{probe}=800`$ $`nm`$$`\mathrm{\Delta }R^{}`$ decreases slowly at low temperature and drops toward zero abruptly around $`T_N`$ (Fig. 2(b)).

The fact, that in both manganites $`\mathrm{\Delta }R^{}`$ vanishes above the transition temperature, is a strong indication that $`\mathrm{\Delta }R^{}`$ arises: a) from quasiparticle excitation involving localized (intragap) states, and b) that the decay of these metastable states involves a spin-flip process. The lifetime $`\tau _{SL}`$ 10 $`\mu `$s of these metastable states is comparable with the spin-lattice relaxation time measured by the $`\mu `$SR technique .

The solid line in Fig. 2(b) shows that $`\mathrm{\Delta }R^{}`$ in LMO follows approximately a $`(T_NT)^\alpha `$ dependence ($`\alpha `$ $``$ 0.5), which is similar to the temperature dependence of the magnetization . As the temperature approaches $`T_N`$ and the magnetic order decreases in LMO the lifetime of the metastable state decreases due to enhanced spin scattering and $`\mathrm{\Delta }R^{}`$ starts to decrease. The non-zero $`\mathrm{\Delta }R^{}`$ above $`T_N`$ indicates that a residual anti-ferromagnetic (AFM) order exists in LMO above the Neel temperature . In contrast to LMO, the $`T`$-dependence of $`\mathrm{\Delta }R^{}`$ is more complex in LCMO. The transient reflectivity first increases with increasing temperature and then for $`T>0.9T_C`$ follows a $`(T_CT)^\beta `$ dependence with $`\beta `$ 0.7 - 0.9 (solid line in Fig. 2(a)). A similar power-law dependence of $`\mathrm{\Delta }R^{}`$ has been observed in other CMR manganites, e.g. Nd0.67Sr0.33MnO3 .

A simple physical model accounts qualitatively for these observations. The various processes giving rise to the photoinduced reflectivity change in LCMO are depicted in Fig. 3 (inset). An ultrashort laser pulse first excites electrons via interband transitions. In our experiments, the films absorbed 1018 \- 1020 photons/cm3 per pulse; comparable to the charge-carrier density ($``$1020-10$`^{21\text{ }}`$holes/cm3) in LCMO; hence, one expects significant electron excitation during ultrashort pump pulse illumination. These hot electrons very rapidly release their energy via electron-electron and electron-phonon collisions reaching QP states near the Fermi energy (step 1). This generates spin waves (magnons) . Spin waves lead to the excitation of metastable states with $`\mu `$s time scale . The QPs can recombine by interaction with states in the pseudogap (step 2) or relax to metastable states via spin-flip processes (magnons are released) caused by strong electron-lattice coupling (step 3). The carriers in the metastable states will relax with a recombination rate $`\gamma =1/\tau _{SL}`$, while magnons are absorbed (step 4). As the pseudogap opens up with increasing temperature and the spectral weight of states at the Fermi level decreases the decay rate of excited QPs decreases (step 2) and more quasiparticles will be scattered into metastable states (step 3). This effect leads to the initial rise of $`\mathrm{\Delta }R^{}`$ which is inverse proportional to the density of states in the pseudogap, i.e., for $`T0.9T_C`$, $`\mathrm{\Delta }R^{}(T_CT)^b`$ with $`b`$ $``$ 0.3 - 0.5 (solid line in Fig. 2(a)). With further increasing temperature ($`T>0.9T_C`$), $`\mathrm{\Delta }R^{}`$ drops significantly toward zero due to the increase of the density of states and occupation of the down-spin conduction electrons, which introduces magnetic disorder.

Figure 4 shows magnetic-field dependent measurements of $`\mathrm{\Delta }R`$ at 305 K (well above $`T_C`$) as a function of time for LCMO thin film grown on NdGaO3 (110) substrate. Magnetic fields are applied in the Faraday geometry, generating an isothermal magnetic entropy change. The results of Fig. 4 are consistent with our model. The quasiparticle relaxation dynamics shows strong magnetic field dependence for $`B`$ ranging from 0.5 to 3 T. The magnetic-field dependence of $`\mathrm{\Delta }R`$ for LCMO is very similar to the temperature dependence of the zero-field $`\mathrm{\Delta }R`$ shown in Fig. 1(a). The increase in the amplitude of $`\mathrm{\Delta }R^{}`$ is due to the increase of $`T_C`$ with applied magnetic filed; this results in an increase of magnetic correlations (spatially and temporally) with enhanced spin alignments. The magnetic correlations are strong enough to overcome thermal fluctuations at the new transition temperature.

The model also correctly explains the temperature dependence of the polaron relaxation time $`\tau _B`$ for LMO and LCMO single crystals (Fig. 3). In LMO single crystal, $`\tau _B`$ is found to be completely temperature independent below $`T_N`$ consistent with a static Jahn-Teller gap. For LCMO single crystal, $`\tau _B`$ remains nearly constant below 0.45 $`T_C`$ and starts to increase above 0.7 $`T_C`$. A similar $`T`$-dependence has been reported in LCMO thin films . The polaron relaxation time $`\tau _B`$ follows a $`(T_CT)^\beta `$ dependence, i.e., $`\tau _B1/\mathrm{\Delta }R^{}`$, shown as a solid line in Fig. 3. This result strongly supports our model that the quasiparticle relaxation dynamics involves fully spin-aligned pseudogap states near $`E_F`$.

In summary, we investigated the spin and quasiparticle relaxation dynamics in LCMO and LMO single crystals and thin films as a function of temperature and magnetic field by time-resolved pump-probe spectroscopy. The quasiparticle relaxation rate is governed by the temperature- and magnetic field- dependent pseudogap in LCMO and by the temperature-independent Jahn-Teller distortion gap in LMO. Both manganites exhibit metastable (localized) states with lifetime $`\tau _{SL}`$ 10 $`\mu `$s, which decay via spin-lattice coupling. Our results show that the coupled dynamics of charge, spin and lattice is strongly correlated with the distinct gap structures in these manganites.

This work is supported in part by NSF through grants: DMR-0137322, IMR-0114124 (WM), and DMR-9876266 (PSU), and the Petroleum Research Fund (PSU).

Current address for C. S. Hong, Department of Chemistry, Korea university, Anam-dong, Sungbuk-ku, Seoul 136-701, Republic of Korea

Figure Captions:

Fig. 1. Transient reflectivity $`\mathrm{\Delta }R`$ of: a). LCMO and b). LMO single crystals around $`T_C`$ (or $`T_N`$) at 800 nm. The dotted lines indicate the zero line.

Fig. 2. Temperature dependence of $`\mathrm{\Delta }R`$ (normalized) at $`\mathrm{\Delta }t=500`$ ps: a) LCMO single crystal ($`\lambda _{pump}=800`$ $`nm/\lambda _{probe}=800`$ $`nm`$) and 400-nm thin film ($`\lambda _{pump}=800`$ $`nm/\lambda _{probe}=5.2`$ $`\mu m`$), and b) LMO single crystal ($`\lambda _{pump}=800`$ $`nm/\lambda _{probe}=800`$ $`nm`$). The solid lines indicate the power-law dependence.

Fig. 3. The relaxation time $`\tau _B`$ as a function of temperature $`T`$ for LCMO and LMO single crystals. The solid line indicates the power-law dependence. The inset depicts a schematic diagram of carrier excitation/relaxation processes in LCMO. For details see text.

Fig. 4. Transient reflectivity $`\mathrm{\Delta }R`$ from LCMO as a function of time for applied fields ranging from 0.5 to 3 T at 305 K ($`\lambda _{pump}`$, $`\lambda _{probe}=800`$ $`nm`$). The dotted lines indicate the zero line.