text:
# Tests for Substructure in Gravitational Lenses

## 1 Introduction

In Dalal & Kochanek (dk02 (2002), DK02 hereafter) we demonstrated that the incidence of anomalous flux ratios in gravitational lenses was consistent with the expected mass fraction of satellites found in cold dark matter (CDM) halo simulations (e.g. Kauffmann Kauffmann93 (1993), Moore et al. Moore99 (1999), Klypin et al. Klypin99 (1999), Bode, Ostriker & Turok Bode01 (2001), Zentner & Bullock Zentner02 (2002)). Our quantitative estimate for the substructure mass fraction was part of an outburst of interest in substructure and anomalous flux ratios of specific lenses (e.g. Mao & Schneider Mao98 (1998), Keeton Keeton01 (2001), Bradac et al. Bradac02 (2002), Chiba Chiba02 (2002), Metcalf Metcalf02 (2002)), in CDM models (e.g. Metcalf & Madau Metcalf01 (2001), Chiba Chiba02 (2002), Zentner & Bullock Zentner02 (2002)), and the difficulty in explaining them based on the properties of the primary lens galaxy (e.g. Metcalf & Zhao Metcalf01 (2001), Keeton, Gaudi & Petters Keetonetal02 (2002), Evans & Witt Evans02 (2002), Möeller, Hewett & Blain Moeller02 (2002), Quadri, Möeller & Natarajan Quadri02 (2002)). In DK02 we argued that alternate explanations to substructure for the anomalous flux ratios were unlikely, but we did not explore the alternate possibilities in detail or develop tests to distinguish substructure from the alternatives. We do so in the current paper.

The alternatives to substructure fall into three broad categories: propagation effects in the interstellar medium (ISM) of the lens galaxy, problems in the “macro” models for the gravitational potential of the lens galaxy, and confusing “microlensing” by the stars in the lens galaxy with the effects of more massive satellites. In this paper we will consider all three of these possibilities. In §2 we study the wavelength dependence of the flux anomalies and develop a simple statistical test to demonstrate that the anomalous flux ratio problem must be due to gravity. In §3 we show that deviations from elliptical shapes appear unable to account for flux anomalies, and that errors in the model for the potential of the primary lens have statistical properties differing from those created by substructure and observed in the data. In §4 we show that it is difficult to explain the anomalous flux ratios of the radio lenses using microlensing. We summarize our results and outline further tests in §5.

## 2 Ruling out the ISM

The ISM can affect flux ratios either through absorption or scattering. In optical/IR observations, many gravitational lenses show wavelength-dependent flux ratios consistent with the effects of dust extinction (e.g. Falco et al. Falco99 (1998)). In the radio, some lenses show the effects of scatter broadening by electrons at low frequencies (e.g. B1933+503, Marlow et al. 1999a ), although in none of these systems is there evidence for a net change in the radio flux. In this section we examine whether the radio lenses show any frequency dependent changes in their flux ratios that could be a signature for the ISM modifying the observed fluxes, and develop a simple, non-parametric test to distinguish the effects of the ISM from the effects of gravity.

### 2.1 The Frequency Dependence of Anomalous Flux Ratios

Almost all mechanisms by which the ISM can modify flux ratios should show frequency-dependent effects, since most scattering processes relevant to the fluxes of radio images become weaker at higher frequencies. For example, weak scintillation causes flux perturbations with a scaling (assuming Kolmogorov turbulence) $`\nu ^{17/12}`$ (Narayan narayan92 (1992)) and the optical depth for free-free absorption scales roughly as $`\nu ^{2.1}`$ (e.g. Mezger & Henderson Mezger67 (1967)). Thus, one way to constrain the ISM as an explanation of the anomalous flux ratios is to examine the allowed frequency dependence of the effect. If the source has an intrinsic spectrum $`f_{\nu ,s}`$ as a function of frequency $`\nu `$, then in the absence of any perturbations from the ISM, the images have spectrum $`f_{\nu ,i}=|M_i|f_{\nu ,s}`$ given the (signed) image magnification $`M_i`$. If the ISM modifies the fluxes through a frequency-dependent optical depth $`\tau _{\nu ,i}`$, then the image fluxes become $`f_{\nu ,i}=|M_i|f_{\nu ,s}\mathrm{exp}(\tau _{\nu ,i})`$.

We examined the optical depth function using the four (or more) image radio lenses with fluxes measured at 5 GHz and either 8 GHz or 15 GHz (MG0414+0534, Katz, Moore & Hewitt Katz97 (1997); B0712+472, Jackson et al. Jackson98 (1998); B1359+154, Myers et al. Myers99 (1999); B1422+231, Patnaik & Narasimha Patnaik01 (2001); B1555+375, Marlow et al. 1999b ; B1608+656, Fassnacht private communication; B1933+503, Sykes et al. Sykes98 (1998); and B2045+265, Fassnacht et al. Fassnacht99 (1999)). We used the same models and magnification estimates as were used to estimate the substructure fraction with a 5% lower bound on the flux errors. The simplest way to illustrate the problem is to assume that one image is strongly affected by the ISM while the others suffer from little or no absorption or scattering. For each image in a lens, we used the fluxes of the other three to estimate the intrinsic spectrum of the source. From this estimate of the intrinsic spectrum we computed the optical depth needed to reproduce the observed flux of the remaining image by computing the necessary optical depth at 5 GHz, $`\tau _5`$, and its spectral index, $`\alpha `$, under the assumption that $`\tau \nu ^\alpha `$ between 5 GHz and the available higher frequency. We estimated the errors solely from the published measurement errors, so they include no uncertainties arising from the lens model or any time variability between the measurement epochs for the different frequencies.

The results are shown in Fig. 1. We have coded the images by their parities (minima versus saddle points) and relative fluxes (brightest versus faintest for each parity) because we will show below that the distributions of the model flux residuals depend strongly on these image identifications. This is seen in Fig. 1 as the concentration of the bright saddle-point images in the direction of higher optical depth compared to the other images. We have dropped four (of 32) images where the estimated optical depth changes sign between the higher and lower frequencies and it is impossible to define a spectral index. Except at zero optical depth, where the slope estimates become unstable and have larger uncertainties, the lenses broadly require $`\alpha 0`$. Since most radio propagation effects have strong wavelength scalings, the ISM seems an unlikely explanation.

### 2.2 A Statistical Test for Substructure

While the spectral indices required to explain the anomalous flux ratios are peculiar for standard scattering or absorption processes, we would prefer to have a test which can distinguish any ISM effect from the gravitational effects of substructure. One approach is to note that the physical properties of the ISM are a local property of the lens, and should have no knowledge of the global properties of the lens geometry. Radial gradients in the properties of the ISM should matter little because in each lens the images of interest effectively lie at the same radius and there is no simple angular correlation between the positions of the images and the major axis of the galaxy. Moreover, creating the anomalous flux ratios depends on the presence of clumped ISM components rather than smoothly distributed components both to create physical conditions extreme enough to have an effect and to differentially affect the lensed images. The only exception to this rule arises for effects such as scintillation or scatter broadening where the effect diminishes as the source size becomes larger – for these cases, the ISM should preferentially perturb the least magnified images because they are the most compact.

The effects of gravity, however, are not determined purely by local conditions because the image positions and fluxes are determined by the competition between the local gravitational potential and the geometric time delay. The magnification tensor, whose eigenvalues determine the image parities, depends on the projected surface density, the projected tidal shear in the gravity, and the redshifts of the lens and the source. As a result, we can attach a global identification to each image based on its parity (maximum, minimum, saddle point of the time delay surface) as well as its magnification with which the ISM should show little (magnification) or no (parity) correlation. Image parities and magnification orderings (but not the precise magnifications) are generic predictions of successful lens models. In particular, for a four image lens, the images alternate parities in the sequence minimum–saddle point–minimum–saddle point as we go around the lens. We can also identify the most and least magnified image of the two minima or saddle points from the geometry of the lens. Unlike the ISM, the effects of substructure depend on the image parity and magnification. The fluxes of highly magnified images are more unstable to small gravitational perturbations than the fluxes of less magnified images (as originally emphasized by Mao & Schneider Mao98 (1998)). Most importantly, the magnification perturbations created by low optical depth substructure have a very strong dependence on the image parity because the perturbations to the fluxes of the saddle points are skewed in the direction of demagnification (Schechter & Wambsganss sw02 (2002), Keeton Keeton02 (2002)). This leads to a simple non-parametric test for distinguishing the effects of the ISM from the effects of gravity – if the flux anomalies depend on parity and magnification they cannot be due to the ISM.

We first demonstrate the effects of substructure on the distribution of flux residuals. We took the same sample of 7 lenses we used in DK02 (MG0414+0534, Hewitt et al. Hewitt92 (1992); B0712+472, Jackson et al. Jackson98 (1998); PG1115+080, Weymann et al. Weymann80 (1980); B1422+231, Patnaik et al. Patnaik92 (1992); B1608+656, Fassnacht et al. Fassnacht96 (1996); B1933+503, Sykes et al. Sykes98 (1998); and B2045+265, Fassnacht et al. Fassnacht99 (1999)). The data were fit using standard macro models (singular isothermal ellipsoids (SIE) plus external shear) to generate an initial model. We then added a substructure mass fraction $`f_{sat}`$ near the images, modeling the substructure as tidally truncated isothermal spheres with critical radius $`b=0\stackrel{}{\mathrm{.}}001`$ (corresponding to approximately $`10^6M_{}`$). With the substructure added we reidentified the lensed images and their fluxes and added astrometry ($`0\stackrel{}{\mathrm{.}}003`$ rms) and flux (10% rms) measurement errors to generate a model data set with substructure. This model data was then refit using the standard macro models, including no weight on the fit to the fluxes of the images. We then compare the observed image fluxes, $`f_{obs,i}`$ to the model image fluxes, $`f_{mod,i}=|M_i|f_{src}`$, found given the magnification $`M_i`$ predicted by the model fitted to the perturbed data. We generated 10 such realizations for each of the 7 lenses. Figure 2 shows the distributions of the residuals, $`\mathrm{log}(f_{obs,i}/f_{mod,i})`$, for the individual image types (saddle point/minimum, brightest/faintest) as a function of the substructure mass fraction, $`f_{sat}`$. As expected, the residual distribution for the brightest saddle point is markedly different from that for all other images unless the substructure fraction is so low as to be masked by the random flux errors.

We now repeat the test for the real data. The same standard macro models were used for each system, and we again assigned no weight to fitting the observed image fluxes. For radio lenses we used the highest frequency measurements from either the VLA or Merlin, and for the optical lenses we used the flux ratios at the longest available wavelength (generally H-band where none of the systems shows significant extinction). This gave us a sample of 8 radio systems (adding B0128+437, Phillips et al. Phillips00 (2000), and B1555+375, Marlow et al. 1999b  to the DK02 sample), 10 optical systems 111The systems are the four-image quasar lenses HE0230–2130 (Wisotzki et al. Wisotzki99 (1999)), HE0435–1223 (Wisotzki et al. Wisotzki02 (2002)), HS0810+2554 (Reimers et al. Reimers02 (2002)), RXJ0911+0551 (Bade et al. Bade97 (1997)), PG1115+080 (Weymann et al. Weymann80 (1980)), H1413+117 (Magain et al. Magain88 (1988)), Q2237+030 (Huchra et al. Huchra85 (1985)), B1422+231 (Patnaik et al. Patnaik92 (1992)), B2045+265 (Fassnacht et al. Fassnacht99 (1999)) and MG0414+0534 (Hewitt et al. Hewitt92 (1992)). The latter 3 systems are also in the radio sample. We used the Castles H-band fluxes if available, otherwise the I-band fluxes.  and 15 joint systems. For our standard model we estimated the source flux using constant fractional errors on the image fluxes. Figure 3 shows the distribution for the residuals, $`\mathrm{log}(f_{obs}/f_{mod})`$ for all presently available quads containing either compact radio sources or bright quasars. Several systems appear in both categories, and for these systems we use the radio fluxes for the joint sample. The distributions of the residuals differ for minima and saddles in the sense that saddle points tend to be fainter than expected and minima tend to be brighter – just as in the simulations of the distributions expected for substructure in Fig. 2.

We used the Kolmogorov-Smirnov (KS) test to estimate the significance of the differences. If we simply compare the distributions of $`\mathrm{log}(f_{obs}/f_{mod})`$ for minima and saddle points, we find that the KS test probabilities for the two distributions to be the same are 2.3%, 28%, 5.5% for the radio, optical and joint samples respectively. Theoretically, we expect the most highly magnified saddle point to have the most discrepant residual distribution because the high magnification makes the flux unstable to perturbations and because the theoretical studies of low optical depth substructure predict that it should show the largest differences (as in Fig. 2). The KS test probabilities for the most magnified saddle point to have the same residual distribution as the other three images is only 0.04%, 5% and 0.3% for the radio, optical and joint samples. In each case, the most highly magnified minimum shows the next lowest probability for a residual distribution agreeing with the other images (1.9%, 59% and 18% respectively), but the significance is lower. In the case of substructure, the magnification perturbation distribution for the minima is harder to differentiate from the measurement uncertainties which probably dominated the residuals for the faintest images.

We can check these significance estimates by examining trials in which the identifications of the images are randomized. We made $`10^4`$ trials in which we assigned the image identifications (1=brightest minimum, 2=faintest minimum, 3=brightest saddle point, 4=faintest saddle point) to the images at random with the restriction that each lens always have the 4 possible image types. We then estimated the probability that these random assignments could produce results similar to that found for the true assignments. The fraction of the time the observed values of the KS statistic were exceeded tracked the KS test probability estimates well. For example, in only 2.3%, 18% and 4.8% of the trials did the KS statistic for the comparison between the distributions of minima and saddle points exceed that of the true distribution. We also examined the distribution of the KS statistic if we bootstrap resampled the lens sample. The KS statistics and probabilities were typical of any sample of lenses drawn (with replacement) from the available 4-image lenses, so the results are not dependent on any particular system or small subset of systems.

We also examined weighting schemes for estimating the source flux other than assuming constant fractional errors for the observed image fluxes. Note, however, that our test is independent of any simple rescaling of the source flux. The KS test makes a non-parametric comparison of the $`\mathrm{log}(f_{obs}/f_{mod})`$ distributions, and a uniform change in the source flux simply shifts the distribution. If we assume constant flux errors rather than constant fractional errors, the probability the saddles and minima have the same distribution is again 2.3%. If we set the source flux using the fluxes of the faintest minimum and saddle, then the probability that the brightest minimum and saddle have the same residual distribution for the radio lenses is 0.02%.

Finally, we can compare the observed distribution to either those expected for Gaussian errors or those we calculated for our Monte Carlo simulations of substructure. For example, while the overall residual distributions of the images are well fit by a log-normal distribution of width $`\sigma _{\mathrm{ln}f}=0.24`$, the residual distributions of the images separately are grossly inconsistent with log-normal distributions. The KS test likelihood for the joint distribution of all images matching a log-normal distribution is $`0.26`$, while the likelihood for the distributions of the individual types matching a single log-normal distribution is $`1.5\times 10^5`$. This is simply another way of quantifying the fact that the residual distributions depend on the image parities and flux orderings. If we compare the residual distributions by image for the radio lens sample to the distributions predicted in the Monte Carlo simulations, we find KS test probabilities for the distributions to be the same of $`4\times 10^7`$, $`2\times 10^5`$, $`2\times 10^4`$, $`0.02`$, $`0.3`$ and $`0.05`$ for the cases with satellite mass fractions of $`f_{sat}=0.002`$, $`0.01`$, $`0.02`$, $`0.05`$, $`0.10`$ and $`0.20`$ respectively. While we do not, at present, feel this is as optimal a way of estimating $`f_{sat}`$ as the methods we used in Dalal & Kochanek (dk02 (2002)), it emphasizes the consistency of the observed dependence of the residual distribution on parity and flux with that expected for substructure.

## 3 Ruling Out Systematic Problems in the Model Potentials

If the anomalous flux ratios cannot be explained by propagation effects in the ISM, the only alternative to some form of substructure in the gravitational field is a systematic problem in the assumed macro model for the lens galaxy. Previous studies of the effects of changes in the macro model show that changes in the radial structure from the standard singular isothermal ellipsoids (SIE) with external shear do not provide a solution (Metcalf and Zhao mz02 (2002), Keeton et al. Keetonetal02 (2002)). This makes physical sense because the images in most of the systems with anomalous flux ratios all lie close to the Einstein radius of the lens, making it difficult for changes in the radial profile to produce changes in the flux ratios. We focus on changes in the angular structure because they can more easily shift the positions of the tangential critical line so as to change the image fluxes.

### 3.1 Lens Models With Deviations from Ellipsoidal Structure

We modified the gravitational potential by adding higher order multipoles of the form

$$\varphi _m=\frac{ϵ_m}{m}r\mathrm{cos}m(\theta \theta _m)$$
(1)
to the primary lens galaxy. These correspond to a surface density

$$\kappa _m=\frac{ϵ_m}{r}\frac{1m^2}{m}\mathrm{cos}m(\theta \theta _m)$$
(2)
and match the expansion used by Evans & Witt (Evans02 (2002)). The new terms have been included in a revised version of the lensmodel package (Keeton Keeton01b (2001)).222The lensmodel package now allows terms of the form $`(ϵ_m/m)r^n\mathrm{cos}m(\theta \theta _m)`$ where the normalization was chosen to reduce to an external shear for $`n=2`$ and $`m=2`$. In lensmodel, the surface density of the SIE model is defined by

$$\kappa =\frac{1}{\sqrt{2}}\frac{b}{r}\frac{q}{\left[1+q^2(1q^2)\mathrm{cos}2(\theta \theta _2)\right]^{1/2}}$$
(3)
where $`q=1ϵ_2`$ is the axis ratio and $`\theta _2`$ and is orientation of the major axis. The critical line, which is an isodensity contour with $`\kappa =1/2`$ for the SIE, has a major axis of $`b`$ and a minor axis of $`bq`$. Thus, if we add an $`m=4`$ term to an SIE model, these normalizations imply a standard amplitude for the deviation of an isophote from the ellipse of $`a_m=ϵ_m(1m^2)/mb`$. Positive $`ϵ_4`$ corresponds to positive $`a_4`$ and “disky” density contours when the $`m=4`$ component is aligned with the quadrupole ($`\theta _4=\theta _2`$). Galaxies with “boxy” isophotes will have negative $`a_4`$ and $`ϵ_4`$.

We will restrict our attention to adding $`m=3`$ and $`m=4`$ terms to the potential of the primary lens galaxy. This should be sufficient to indicate their potential importance without the models being too under-constrained. Observational evidence indicates that the amplitude of these higher order multipoles is small. For example, Bender et al. (bender (1989)) and Rest et al. (Rest01 (2001)) examine a sample of local ellipticals and quantify the deviations of the isophotes from ellipses. They find that the deviations are dominated by the $`m=4`$ term, $`a_4\mathrm{cos}(\theta \theta _4)`$, with a typical amplitude of $`|a_4|0.01`$. Deviations from ellipticity in projections of simulated halos are of comparable amplitude (Heyl et al. heyl (1994), Burkert & Naab Burkert03 (2003)). Considerably less is available on the distribution of misalignment angles, $`\theta _m\theta _2`$, between the higher order multipoles and the dominant ellipse. We know from the statistics of lens models that the major axes of the mass and the light are closely aligned (Kochanek Kochanek01b (2002)). The effect of these higher order terms on lenses were first studied by Evans & Witt (ew (2001)), who pointed out that larger values of $`a_4`$ than are typically observed would lead to significant cross sections for lenses producing more than 4 observable images.

We fit the 7 lenses in the DK02 sample with a sequence of four models whose results are presented in Table 1. In the first model (“astrometry”) we fit only the image positions using our standard SIE plus external shear model. In the second model (“astrometry$`+`$flux”) we added the constraints from the image fluxes assuming 5% uncertainties. In either case, 6 of the 7 lenses have unacceptably high $`\chi _{flx}^2`$ fit statistics for the image fluxes even though the standard models provide good fits to the astrometric data, $`\chi _{tot}^2\chi _{flx}^2`$. Only B1608+656 has a $`\chi _{flx}^2`$ largely consistent with 5% measurement uncertainties. The failure of these models to reproduce the flux ratios is the origin of the entire problem of flux ratio anomalies. In the third model we added the $`m=3`$ and $`m=4`$ multipoles to the potential of the primary lens, but constrained them to be aligned with the ellipsoid ($`\theta _m=\theta _2`$ but allowing any sign for $`ϵ_m`$). For MG0414+0534, B0712+472, PG1115+080 and B1933+503 the new terms lead to a significant reduction in $`\chi _{flx}^2`$, but in no case does it become a statistically acceptable fit. The amplitudes of the multipoles are unreasonably large for B0712+472 and B1933+503. Finally, in the fourth model, we allowed the $`m=3`$ and $`m=4`$ multipoles to be misaligned with respect to the ellipsoid. In these models the flux ratio anomalies implied by the standard lens models for PG1115+080, B1422+231 and B1933+503 are gone, but the multipole amplitudes are still unreasonably large.

Thus, like Evans & Witt (Evans02 (2002)), we find that adding higher order multipoles to the potential could explain some of the flux ratio anomalies (3 of the 6 in the DK02 sample), if the positions and fluxes of the 4 compact images were the only available data. The surface density contours implied by these models are not very physical, as shown for PG1115+080 and B1933+503 in Figs. 5 and 5. The improved fit with the higher order multipoles is easily understood for B1933+503. In this lens, image 4 is anomalously bright for an ellipsoidal model, and the new models solve the problem by making the tangential critical line lie much closer to image 4 so that it can be more highly magnified. The cause of the improvement in PG1115+080, where the anomaly is due to the flux ratio between images A1 and A2, is less obvious. In other cases, like MG0414+0534, the surface density contours, while very boxy, look more reasonable (Fig. 6).

For three of the DK02 lenses with anomalies there are additional lensed components beyond the four compact images which can be used to test these models. In MG0414+0534 there are VLBI jets associated with each of the quasar images (Ros et al. Ros00 (2000), Trotter, Winn & Hewitt Trotter00 (2000)), in PG1115+080 there is the Einstein ring image of the quasar host galaxy (Impey et al. impey98 (1998), Kochanek et al. Kochanek01 (2001)), and in B1933+503 there are additional multiply-imaged components of the radio source (Sykes et al. Sykes98 (1998), Cohn et al. Cohn01 (2001), Munoz, Kochanek & Keeton munoz01 (2001)). Table 2 shows the results of fitting the same series of models to these three lenses with the addition of the new constraints. For MG0414+0534, where the models were already unable to explain the anomalous flux ratios, little changes with the addition of the VLBI constraints, although the amplitudes of the multipoles change as shown in Fig. 6. For the two lenses where it appeared possible for the higher order multipoles to eliminate the anomalies, we find that the new constraints eliminate these solutions. Figs. 5 and 5 show the surface density contours for these new models, and it is immediately clear that the added constraints force the lenses to be significantly more regular and ellipsoidal.

### 3.2 Statistical Checks of the Macro Model

Thus, we find that higher order multipoles can explain only 1 (B1422+231) of the 6 anomalous flux ratio lenses in the DK02 sample given the available data. Keeton et al. (Keetonetal02 (2002)) also noted that B1422+231 was less problematic than many of the other systems. Nonetheless, the perturbation amplitudes needed to explain B1422+231 are much larger than observed in real galaxies, halo simulations or in the lenses where additional constraints allow us to determine $`a_3`$ and $`a_4`$. We conclude that higher order multipoles are not a likely explanation, but we would like a more generic test of this conclusion. The basic problem is that most lenses provide only enough data to constrain the standard ellipsoidal models, so more complex models that can explain the image fluxes are only constrained by the degree to which the models are viewed as physical.

We noted in §2.2 that there is a statistical tendency for the anomalous flux ratios to appear as a suppression of the flux of the brightest saddle points relative to the predictions of smooth models, and that this was a characteristic property of low optical depth substructure. It is likely that any specific model which can explain the anomalies in individual lenses will be unable to reproduce this statistical property for randomly selected lenses. The worst case for examining this expectation is the effect of the $`m=4`$ multipole, because it is the largest amplitude higher order multipole with a symmetry matching a common image configuration, the cruciform quads like Q2237+0305. If an $`m=4`$ multipole is aligned with the ellipsoid it will preferentially affect the saddles as compared to the minima of a cruciform quad, but it will symmetrically affect both saddles and both minima. If it is randomly misaligned with respect to the ellipsoid, then it will not distinguish between saddles and minima, and if the observed lens geometries are dominated by merging pairs (fold catastrophes, like PG1115+080) or merging triples (cusp catastrophes, like B1422+231), then it will again be unable to distinguish saddles from minima.

We tested these hypotheses for models consisting of an ellipsoid with an aligned (boxy or disky) $`a_4`$ component or a randomly oriented $`a_4`$ component for a range of models for the lens population or typical lens geometry. If we simulate $`m=4`$ models with the amplitudes typical of real elliptical galaxies or simulated CDM halos ($`|a_4|0.02`$), the resulting $`5\%`$ flux perturbations are too small to be relevant for explaining the anomalous flux ratios. We must use models of abnormally large amplitude, $`|a_4|=0.05`$, to produce sufficiently large perturbations. After generating an image configuration, we fit the images using our standard SIE $`+`$ shear models, fitting only the astrometric data and ignoring the fluxes. We then examined the statistics of the flux residuals after dividing the images into magnification and parity subsamples just as in §2.2. We focused on cusp and fold configurations defined by lenses with an opening angle between two images of $`<30^{}`$. Of the DK02 sample, the systems MG0414+0534, B0712+472, PG1115+080, B1422+231, and B2045+265 would be classified as fold or cusp under this definition, while B1608+656 and B1933+503 just fail to make the cut.

The results of our Monte Carlo calculation are displayed in Figure 7a-c. As expected, $`a_4`$ perturbations to fold or cusp lenses do not systematically distinguish between image parities for disky, boxy, or randomly oriented $`m=4`$ multipoles. This is in stark contrast to the properties of real fold and cusp lenses where the bright saddle points are preferentially demagnified compared to the other images (Fig. 7d). Clearly, $`m=4`$ perturbations are incapable of producing this behavior, with KS test probabilities of less than $`0.002`$ than any of these distributions agree with the distribution observed for the cusp/fold lenses in DK02. We experimented with still higher order multipoles ($`m>4`$), and found that they also fare poorly in both systematically demagnifying saddles and differentiating between saddles and minima. The differences between distributions found in the Monte Carlo simulations and in the real data mean that the flux ratio anomalies cannot be due to higher order multipoles. This would be true even if the models with higher order multipoles could successfully model the anomalous flux ratios in every lens!

## 4 Ruling Out Microlensing of Radio Sources

We are left only with substructure in the gravitational potential as a viable explanation of the phenomenon. Real lenses have substructures on two mass scales, stars and satellite halos, so the last step of our argument is to rule out stars as the source of the anomalous flux ratios of lensed radio sources. Microlensing is certainly an important phenomenon in the optical, where its effects are directly observed as uncorrelated time variability in lensed quasars (e.g. Q2237+0305, Wozniak et al. Wozniak00 (2000)). The principal arguments against significant microlensing in radio lenses are laid out by Koopmans & de Bruyn (koopmans (2000)). The key point is that the typical scale of the caustic network is set by the Einstein radii of the stars ($`\theta _\mathrm{E}=3[M/M_{}]^{1/2}\mu `$as for a star of mass $`M`$, a lens at $`z_l=0.3`$ and a source at $`z_s=2`$, hence the name microlensing), while the Compton limit for the minimum possible radio source size is given by

$$S(\nu )=\frac{2kT_b\nu ^2}{(1+z)c^2}\pi \mathrm{\Delta }\theta _{\mathrm{min}}^2.$$
(4)
Table 3 shows that the minimum source sizes for the four-image CLASS lenses are typically $`\theta _{\mathrm{min}}>10\mu `$as, or roughly 10 times larger than $`\theta _\mathrm{E}`$ for a stellar mass function dominated by low-mass stars. A sufficiently non-thermal electron distribution can change these limits, but it is unlikely to reduce them by an order of magnitude. While microlensing leads to fractional flux variations of order unity for sources smaller than the typical caustic size, larger sources average over the magnification pattern and show far smaller variations. For sources larger than $`20\mu `$as we would expect $`5\%`$ fluctuations from stellar microlensing. For microlensing to produce the $`20\%`$ flux anomalies of the radio lenses we would need to have source sizes of order $`\mu `$as.

One means of escaping the Compton limit on the radio source sizes is to make most of the radio flux come from a relativistically beamed source. If a significant fraction of the source flux comes from a source moving towards us at relativistic speeds, then the brightness temperature of the beamed component can be increased by a Doppler factor $`𝒟=[\gamma (1\beta \mathrm{cos}i)]^1\gamma `$, allowing a corresponding reduction in the source size. For example, there are $`5\%`$ uncorrelated changes in the radio fluxes of the lensed images in B1600+434 that can be explained by microlensing a beamed radio component with $`\gamma 10`$ (Koopmans & de Bruyn koopmans (2000)). Much larger Lorentz factors would be needed for beaming to explain the $`20\%`$ amplitude of the flux anomalies.

In the absence of beaming, microlensing leads to changes in the image flux ratios on time scales of years, while CDM substructure would lead to changes only on time scales of millenia or longer. If there is beaming, so that the characteristic velocities are superluminal rather than simply the peculiar velocities of the lens, source and observer, then the time scales will be $`1000`$ times shorter. Where radio sources have been monitored for extended periods (e.g. B1600+434, Koopmans & de Bruyn koopmans (2000); B1933+503, Biggs et al. biggs (2000); B1608+656, Fassnacht et al. Fassnacht02 (2002)) the average flux ratios show little time dependence. While this is only weak evidence against microlensing in the absence of beaming, it is strong evidence against microlensing if there is beaming because of the very short time scales created by the relativistic motions. Finally, if microlensing were important, we would expect a significant frequency dependence to the flux ratios because we expect the source size to be a function of wavelength. As discussed in §2.1, this frequency dependence has not been observed.

Because of the stress we have put on the dependence of the flux residuals on image parity, we examined whether the differences between saddle points and minima we discussed earlier persists for source sizes comparable to the Einstein radii of the substructure. We ran inverse ray-tracing simulations (e.g. Schneider et al. Schneider92 (1992)) using 102402 rays and the particle-mesh (PM) force calculation algorithm (Hockney and Eastwood hockney (1981)). The PM method enjoys several advantages for the microlensing problem, chiefly in the simplicity and speed of the algorithm. Using this code, we simulated microlensing at a saddle-point and at a minimum, each with an unperturbed macro-model magnification of about 10. The saddle has a total convergence and shear (coming from the macro model) of $`\kappa _0=0.525,\gamma _0=0.575`$, while the minimum has $`\kappa _0=0.475,\gamma _0=0.425`$. We lay down a random distribution of microlenses with Einstein radii $`\theta _\mathrm{E}=1\mu `$as and a number density corresponding to convergence $`\kappa _{}=0.075`$. The remaining surface density, $`\kappa _0\kappa _{}`$, was in a smooth component.

In Figure 8 we plot the distribution of flux residuals as a function of source size. For a source small compared to $`\theta _\mathrm{E}`$, like the case $`r_s=0.13\mu `$as corresponding to our grid resolution, we see that the magnification distribution for the saddle point is skewed towards demagnification compared to that for the minimum, in agreement with previous results (e.g. Schechter & Wambsganss sw02 (2002); see their Fig. 3). For larger sources, which we illustrate for the cases $`r_s=\theta _\mathrm{E}=1\mu `$as and $`3\theta _\mathrm{E}=3\mu `$as, the widths of the distribution narrow rapidly (with $`\sigma _\mu /\mu _01/r_s`$ roughly, e.g. Koopmans & de Bruyn koopmans (2000)) and the differences between the distributions for saddle points and minima vanish. The mean $`\mathrm{log}_{10}\mu /\mu _0`$ provides a measure of the skewness, and for source sizes of $`r_s=0.13`$, $`0.56`$, $`1`$ and $`3\mu `$as we find values of $`\mathrm{log}_{10}\mu /\mu _0=0.26`$, $`0.025`$ $`0.01`$ and $`0.002`$ respectively. Evidently, even source sizes as small as $`\theta _\mathrm{E}/2=0.5\mu `$as are sufficient to erase the asymmetry in the distribution of perturbations to saddle-point images. In reality, radio QSO’s are expected to have much larger angular sizes (10–30$`\mu `$as, see Table 3), giving even narrower, more symmetric distributions than those plotted in Figure 8.

## 5 Summary and Discussion

We have considered three alternatives to substructure as the source of the flux anomalies observed in four-image radio lenses: propagation effects like refractive scattering in the interstellar medium (ISM) of the lens galaxy, problems in the macro model describing the gravitational potential of the lens galaxy, and microlensing by ordinary stars in the lens galaxy. We find that none of these alternatives are likely to explain the phenomenon.

Generating the flux anomalies using the ISM fails on two counts. First, the anomalies in the flux ratios of the radio lenses are essentially independent of wavelength, unlike any normal propagation effect in the ISM (scintillation, refractive scattering, free-free absorption). Any attempt to use the ISM as an explanation of the anomalous flux ratios must find a mechanism to produce a nearly frequency independent optical depth – the radio equivalent of the “gray dust” proposals for Type Ia supernovae (Aguirre Aguirre98 (1998)). Second, the flux ratio anomalies depend on the image parity and magnification. In particular, the brightest saddle point image has a tendency to be demagnified relative to the other images, which is a characteristic of low optical depth substructure (Schechter & Wambsganss sw02 (2002), Keeton Keeton02 (2002)). Since the ISM should preferentially perturb the least magnified images, which will be the most compact, and show no dependence on the image parity, we can statistically rule out any propagation effect.

Solving the flux anomalies by modifying the smooth gravitational potential of the lens also fails on two counts. Following other studies (Evans & Witt Evans02 (2002), Möeller, Hewett & Blain Moeller02 (2002), Quadri, Möeller & Natarajan Quadri02 (2002)), we investigated adding higher order multipoles ($`\varphi \mathrm{cos}m\theta `$ with $`m>2`$) to the potential. For the models we attempted (which were restricted to adding $`m=3`$ and $`m=4`$) we could find solutions without flux anomalies for only one of the six lenses from DK02 which had anomalies in our standard ellipsoidal lens models. The fits required the higher order multipoles to be significantly misaligned with respect to the axes of the ellipsoid and had perturbation amplitudes larger than seen in real galaxies or halo simulations. In three of the lenses, the presence of additional constraints, either Einstein ring images of the quasar host galaxy or additional lensed components of the radio source, allowed us to measure directly the magnitude of the higher multipoles in the gravitational potential. In every instance, the $`m=3`$ and $`4`$ multipoles were constrained to the small amplitudes typically observed in galaxies and halo models, largely justifying the use of ellipsoidal mass models to describe lens galaxies.

The second problem with the higher order multipole models is that they generally cannot reproduce the observed parity dependence of the flux anomalies. They can do so in very specific cases where the symmetry of the lens (a cruciform quad like Q2237+0305) matches the symmetry of the multipole (a $`\mathrm{cos}4\theta `$ multipole), and the multipole has very large amplitude and is aligned with the dominant ellipsoidal potential. However, none of the radio quads with anomalies have this symmetry (they generally have merging image pairs or triples associated with fold and cusp caustics) and the models require multipoles that are not aligned with the ellipsoid. If we examine the distribution of anomalies found by fitting lenses with higher order multipoles using only standard ellipsoidal lens models, we find that the statistical properties of the brightest saddle point do not stand out from those of the other images. This allows us to rule out these models as an explanation even if they could fit the anomalies in each individual lens (which they cannot).

The failure of the ISM or changes in the smooth potential to explain the data means that the explanation must be substructure in the potential. These substructures must either be dark objects, since luminous satellites are not abundant enough to account for the preponderance of anomalous fluxes (e.g. Chiba Chiba02 (2002)), or the stellar populations in the lens galaxy. Either source of substructure will lead to the statistical differences between saddle points and minima if the optical depth is low (Schechter & Wambsganss sw02 (2002), Keeton Keeton02 (2002)). The angular scales of the magnification patterns are very different for stars (microarcseconds) and satellites (milliarcseconds) and by using radio sources we should be averaging out the stellar magnification patterns to leave only the contribution from the satellites. We tested this by taking typical microlensing patterns for saddle point and minimum images, which show the expected differences for point sources, and smoothing on larger and larger scales. As expected, the magnification distributions for the saddle point and the minimum show significant differences only when the source size is smaller than the Einstein radius of the stars. Source sizes of even one microarcsecond are sufficient to eliminate the differences, so radio sources could have angular sizes even an order of magnitude smaller than the Compton limit without microlensing making a significant contribution to the flux anomalies of radio lenses.

In summary, CDM substructure remains the best explanation of the flux ratio anomalies of gravitational lenses. The most powerful piece of evidence is the statistically significant difference between the anomalies in saddle points and minima and between the bright and faint images. Substructure makes a very specific prediction that the brightest saddle point should show a very different distribution from the other images, as observed in the data, that is very difficult or impossible to reproduce using either the ISM or changes in the smooth potential.

However, the issues raised by these considerations point to future observations that can further clarify the origin of anomalous fluxes in lens systems. ISM effects can be more strongly constrained by measuring flux ratios at still higher frequencies (e.g. 43 GHz at the VLA). Mid-infrared (5–10$`\mu `$m) flux ratios, where the wavelength is far too short to be bothered by electrons and far too long to be bothered by dust, are difficult to measure but completely insensitive to the ISM (e.g. Agol, Jones & Blaes Agol00 (2000)). The mid-infrared source sizes should also be large enough to be affected only by satellites, adding a further check for whether we can be misinterpreting microlensing effects as substructure. Integral field spectroscopy, to compare the flux ratios in the optical continuum to those in the broad emission lines, provides another way to separate the effects of microlensing and substructure because the broad line emitting regions should be significantly larger than the source of the optical continuum (e.g. Moustakas & Metcalf Moustakas02 (2002)). Monitoring the lenses, either in the radio or the optical, to look for changes in the flux ratios, can also be used to distinguish microlensing effects from CDM substructure. Observations to find additional lensed structures are the best approach to determining whether more complicated lens potentials are needed. Clean constraints on more complicated angular structures can be obtained by analyzing the shapes of the Einstein ring images of host galaxies found in deep, high-resolution infrared images (see Kochanek et al. Kochanek01 (2001)).

We thank Roger Blandford, Scott Gaudi, Gil Holder, Chuck Keeton, Feryal Özel and David Rusin for helpful discussions. We also thank Chuck Keeton for rapidly adding the higher order multipole potentials to the lensmodel package, and Chris Fassnacht for giving us the unpublished 5 and 15 GHz flux measurements for B1608+656. ND acknowledges the support of NASA through Hubble Fellowship grant #HST-HF-01148.01-A awarded by STScI, which is operated by AURA for NASA, under contract NAS 5-26555. CSK is supported by the Smithsonian Institution and NASA grant NAG5-9265.


format: text

text:
# Optical and NIR monitoring of the GRB020405 afterglow

## 1. Introduction

GRB020405 was detected by the Interplanetary Network on 2002 April 5.02877 UT with a duration of $``$40 s (Hurley et al. 2002). Its optical afterglow was detected by Price et al. (2002a) and confirmed by subsequent observations (e.g. Covino et al. 2002). Optical spectroscopy allowed Masetti et al. (2002a) and Price et al. (2002b) to determine the redshift of the GRB, $`z`$ = 0.691.

Here we report on optical imaging, spectroscopy and polarimetry, and near-infrared (NIR) imaging of the GRB020405 afterglow acquired as part of our ongoing programs of optical/NIR follow-up of GRB afterglows at ESO and other telescopes. We also included the available archival HST pointings. In particular, our NIR observations allowed the first detection of this afterglow at these wavelengths. A more detailed presentation will appear in Masetti et al. (2003, in preparation).

## 2. Observations

Optical and NIR photometry has been accomplished at several ESO telescopes over a period of 10 days, starting on 6 April 2002. Early optical imaging was acquired also at TNG and WHT (Canary Islands, Spain) and at the 1.04-m UPSO (Naini Tal, India). The ground-based data set has been complemented by archival HST data taken at later epochs.

Optical spectra were taken at VLT-Melipal on April 6 and 7 with a dispersion of 5.5 and 2.6 Å/pix, respectively. In the second spectroscopic observation the slit was rotated of about 40 with respect to the North-South direction in order to include in the slit the OT and two nearby galaxies (see Sect. 3). Furthermore, a $`R`$-band polarimetric measurement was acquired at VLT-Melipal on April 6.

## 3. Results

Light curves. In order to take into account the contribution of the host galaxy to the optical transient (OT) emission in the ground-based photometry, we followed two ways: (i) we fitted the light curves with a simple power law plus constant; (ii) we measured, on the latest available HST images, the contribution of the host galaxy complex within an aperture radius matching the ground-based telescopes PSFs. The two methods gave consistent results, and showed that the host contribution to the OT luminosity is of the order of few percent (and thus negligible within the uncertainties) in the first days of observation.

The host-subtracted ground-based optical measurements are fitted by a single power law decay, with index $`\alpha `$ = 1.54$`\pm `$0.06. The host contribution in the HST images was estimated and subtracted as well by applying a PSF-based method for the OT removal. The addition of the HST data points to the $`VRI`$ light curves made it apparent the presence of a deviation from the early single power-law behaviour, i.e. a slight rebrightening followed by a decay (Fig. 1, left). In the $`RIJ`$ light curves, this “bump” is equally well fitted with an emerging SN akin 2002ap, and 1.3 mag brighter, at the OT redshift ($`z`$ = 0.691), and with a SN1998bw dimmed by 0.6 mag (see also Price et al. 2002b). $`V`$-band data (Fig. 1, upper left) are instead poorly fitted by either SN. Alternatively, the bump may be interpreted with a shell collision re-energization scenario, as proposed by Kumar & Piran (2000) and Beloborodov (2002). This explanation would be supported by the satisfactory fit of the light curves past the bump (from day 20 to 70 after GRB) with a second power law of index $`\alpha ^{}`$ = 1.85$`\pm `$0.15 (Fig. 1, right).

Spectra. The spectra acquired on April 6 and 7 with VLT-FORS1 show Balmer, \[OII\] and \[OIII\] lines in emission and FeII and MgII in absortion. From these, we measure a redshift $`z`$ = 0.691$`\pm `$0.002 for the GRB. Both absorption and emission features in the OT spectrum are consistent with this value. To acquire both the spectrum of the OT and those of Galaxies 1 and 2 located at 2′′ and $``$6′′ southwest of it, on April 7 the 1′′-wide FORS1 slit was rotated by 40 towards East with respect to the N-S direction (Fig. 2, left panel).

From the detection of \[OII\] and \[OIII\] emission lines in their spectra, both Galaxies 1 and 2 appear to be at a substantially lower redshift, $`z`$ = 0.472$`\pm `$0.002, than the OT host galaxy (Fig. 2, right). Thus, although angularly close to the host galaxy complex of GRB020405, Galaxy 1 is not interacting with the GRB host galaxy, as formerly proposed by Masetti et al. (2002b). The preliminary estimate of the redshift of Galaxy 1 was incorrect, due to an improper subtraction of the host complex contribution from its spectrum.

Optical-NIR spectral flux distribution. Using photometry data, we have constructed 5 optical-NIR broadband spectra. The data points were corrected for the Galactic absorption ($`E(BV)`$ = 0.055; Schlegel et al. 1998; Cardelli et al. 1989) and converted into fluxes according to Fukugita et al. (1995) for the optical and to Bersanelli et al. (1991) for the NIR. The spectra of the first two epochs were not corrected for the host contribution because this is known for the $`VRIJ`$ bands only; however, given this was quite modest at those epochs, we simply added a 5% error in quadrature to the uncertainties in the optical-NIR fluxes.

The two earliest spectra clearly show a break at 2.5$`\times `$1014 Hz, i.e. in the $`J`$ band, with slopes $`\beta _{\mathrm{NIR}}`$ = 0.65$`\pm `$0.2 and $`\beta _{\mathrm{opt}}`$ = 1.3$`\pm `$0.2. This break can be interpreted as the synchrotron cooling frequency $`\nu _\mathrm{c}`$ in the simplest case of a spherical fireball expanding in a homogeneous medium (Sari et al. 1998) and with electron distribution index $`p`$ 2.6. Broadband spectra of the three following epochs, made with $`VRI`$ points, were instead plotted by subtracting the host contribution. The spectral slopes on these three epochs are consistent with the optical one on the first two epochs. The broad-band spectra of the OT emission during and after the “bump” have a steeper power law shape, $`\beta `$ = 3.5$`\pm `$0.5, in agreement with the findings of Price et al. (2002b).

Polarimetry. Our $`R`$-band polarimetry data indicate for the OT a linear polarization $`P_{\mathrm{OT}}`$ = 1.5$`\pm `$0.4 % and a polarization angle $`\theta _{\mathrm{OT}}`$ = 172$`\pm `$8. This result is corrected for possible instrumental and interstellar polarization by using field stars and polarization standard stars. This value for the polarization is consistent with those measured by Covino et al. (2002), but is at variance with that ($`P`$ 10%) obtained by Bersier et al. (2002) from $`V`$-band observations acquired nearly simultaneously with ours.

## References

Beloborodov, A.M. 2002, ApJ, submitted (astro-ph/0209228)

Bersanelli, M., Bouchet, P., & Falomo, R. 1991, A&A, 252, 854

Bersier, D., et al. 2002, ApJ, submitted (astro-ph/0206465)

Cardelli, J.A., Clayton, G.C. & Mathis, J.S., 1989, ApJ, 345, 245

Covino, S., et al. 2002, A&A in press (astro-ph/0211245)

Fukugita, M., Shimasaku, K., & Ichikawa, T. 1995, PASP, 107, 945

Hurley, K., Cline, T., Frontera, F., et al. 2002, GCN #1329

Kumar, P., & Piran, T. 2000, ApJ, 532, 286

Masetti, N., Palazzi, E., Pian, E., et al. 2002a, GCN #1330

Masetti, N., Palazzi, E., Maiorano, E., et al. 2002b, GCN #1375

Price, P.A., et al. 2002a, GCN #1326

Price, P.A., et al. 2002b, ApJ, submitted (astro-ph/0208008)

Sari, R., Piran, T., & Narayan, R. 1998, ApJ, 497, L17

Schlegel, D.J., Finkbeiner, D.P., & Davis, M. 1998, ApJ, 500, 525


text:
# A composite open resonator for compact X-ray source

## 1 Introduction

Currently there has been considerable progress in the development of super-reflection mirrors, high-finesse optical resonators (super-resonators) and frequency stabilized high power cw and pulsed mode-locked lasers. Optical resonators can be used as a photon storage to accumulate a very high laser power. By this method, a finesse of the optical resonators of $`10^6`$ can be achieved, which means the laser power can be enhanced by this order . The high reactive power stored in the super-resonators can be used in gravitational experiments in astrophysics, Laser-Electron Storage Rings  \- , gravitational analogue of lasers (grasers)  or other applications.

The main reason for mirror degradation of high finesse dielectric open resonators of free-electron lasers based on storage rings is the deposition of chemical components on the mirror surface by X-ray synchrotron radiation . The same degradation can occur in Laser-Electron Storage Rings where the backward Compton scattering of laser photons by electron beams is used for production of X-ray and $`\gamma `$-ray radiation , . Possible solution of the mirror degradation problem is using of a composite resonator. Below the results of calculation of the dependence of the finesse of such a composite resonator on its parameters are presented and analyzed.

## 2 Quality of the composite open resonator

Let a composite open resonator consist of two mirrors (see Fig.1). The reflectivity of the first mirror $`M_1`$ and the main part of the second mirror $`M_2`$ are $`R_1`$. The second mirror has a circular insertion at the axis of the resonator made from another material with the reflectivity $`R_2`$. The radius of the insertion is ”$`a`$”, the rms radius of the fundamental mode of the resonator at the surface of the mirror is ”$`\sigma _{Lm}`$”. Both mirrors can be dielectric multilayer mirrors and the insertion part in one in two mirrors can be coated by $`Au`$.

Fig 1: A composite open resonator.

The intensity of the Gaussian laser beam is distributed by the law

$$I_L=\frac{P_0}{2\pi \sigma _L^2}e^{\frac{r^2}{2\sigma _L^2}},$$
(1)
where $`P_0=I_L𝑑S`$ is the initial power of the laser beam stored in the resonator; $`dS`$, the element of the area; $`\sigma _L=\sigma _{L\mathrm{\hspace{0.17em}0}}\sqrt{1+s^2/l_R^2}`$, the dispersion of the laser beam in a point $`s`$; $`\sigma _{L\mathrm{\hspace{0.17em}0}}=\sigma _L(s=0)`$; the point $`s=0`$ corresponds to the waist of the laser beam; $`l_R=4\pi \sigma _{L\mathrm{\hspace{0.17em}0}}^2/\lambda _L`$, the Rayleigh length; $`\lambda _L`$, the laser wavelength.

In a composite resonator a part $`P_1=P_0R_1`$ of the laser power is reflected by the first mirror. A part $`(P_1\mathrm{\Delta }P_1)R_1`$ of this power is reflected by external part of the second mirror and the another part $`\mathrm{\Delta }P_1R_2`$ by the internal one, where $`\mathrm{\Delta }P_1=P_1[1\mathrm{exp}(a^2/2\sigma _{Lm}^2)]`$ is the energy incident upon the insertion part of the second mirror. The round trip reflected power $`P_2=(P_1\mathrm{\Delta }P_1)R_1+\mathrm{\Delta }P_1R_2`$ can be presented in the form

$$P_2=P_0R_1^2[1\frac{R_1R_2}{R_1}(1e^{a^2/2\sigma _{Lm}^2})].$$
(2)
After $`n`$ round trips the laser light power can be presented in the form

$$P(t)P_0(\frac{P_2}{P_0})^n=P_0e^{n\mathrm{ln}(P_2/P_0)}=$$

$$P_0e^{n\mathrm{ln}\{R_1^2[1\frac{R_1R_2}{R_1}(1e^{a^2/2\sigma _{Lm}^2})]\}},$$
(3)
where $`n=ct/2L`$; $`L`$, is the resonator length.

The finesse of the open resonator can be determined by the equation $`F=(2\pi /T)[P/(P/t)]`$. In this case, according to (3), the finesse of the composite resonator can be presented in the form

$$F=\frac{2\pi }{\mathrm{ln}\{R_1^2[1\frac{R_1R_2}{R_1}(1e^{a^2/2\sigma _{Lm}^2})]\}}=$$

$$F_0\frac{\mathrm{ln}R_1^2}{\mathrm{ln}\{R_1^2[1\frac{R_1R_2}{R_1}(1e^{a^2/2\sigma _{Lm}^2})]\}},$$
(4)
where $`F_0=F|_{R_1=R_2}=\pi /\mathrm{ln}R_1\pi /(1R_1)`$.

According to (4), in the approximations $`(1R_1)/(1R_2)1`$, $`1R_11`$, $`1R_21`$, $`a/\sigma _{Lm}1`$, the finesse of the resonator will be decreased $`F_0/F`$ times when the radius

$$a\sqrt{2}\sigma _{Lm}\sqrt{\frac{1R_1}{1R_2}(\frac{F_0}{F}1)}.$$
(5)
Example. The relativistic factor of electrons $`\gamma =10^3`$, $`\lambda _L=10`$ $`mkm`$, $`\sigma _{L\mathrm{\hspace{0.17em}0}}=50`$ $`mkm`$, $`L=2m`$, $`F/F_0=0.5`$, $`(1R_2)/(1R_1)=10^2`$.

In this case, according to (1), (5), the Rayleigh length $`l_R=3.14`$ $`mm`$, the mode size $`\sigma _{Lm}=15.9`$ $`mm`$, the radius of the X-ray beam at the mirrors $`\sigma _{Xray}L/2\gamma =1`$ $`mm`$, the radius of the insertion $`a=0.1\sqrt{2}\sigma _{Lm}2.25`$ $`mm`$ $`>\sigma _{Xray}`$. The magnification $`M`$, defined as the ratio of the mode size on the mirrors to the mode size at the waist in the resonator, must be chosen to be high ($`10^2÷10^3`$) as a compromise between mirror degradation, resonator finesse and source issues.

## 3 Conclusion

We have considered a composite open resonator. One of its mirrors has a small circular insertion at the axis of the resonator made from another much lower finesse material. It was shown that if the dimension of the insertion is larger than the effective transverse dimension of the X-ray beam and much less than the diameter of the fundamental mode the mirror degradation time can be increased to a great extent and the finesse of the resonator is decreased unessentially. Such conditions can be fulfilled when the magnification $`M1`$.

This work was supported partly by the Russian Foundation for Basic Research, Grant No 02-02-16209.


text:
# The wavelength of neutrino and neutral kaon oscillations

## 1 Introduction

When neutral particles are produced in a flavor eigenstate that is not also a mass eigenstate, the resultant system is, in general, a superposition of mass eigenstates. If so, the system may oscillate between the different flavor eigenstates. This situation has been familiar for many years for the case of neutral kaons; if these are produced in one of the strangeness eigenstates, $`K^0(S=1)`$ or $`\overline{K}^0(S=1)`$, the system is a superposition of the mass eigenstates, $`K_L`$ and $`K_S`$, and oscillates between the two strangeness eigenstates. Recently, strong evidence has been found that neutrinos show a similar behaviour .

The standard quantum-mechanical treatment of kaon oscillations  has been known for many years and results in an expression relating the wavelength of the strangeness oscillations to the mass difference between the mass eigenstates, $`\delta m=m_Lm_S`$. In the last decade, several papers have appeared which question this treatment, sometimes resulting in a different relation between the wavelength and $`\delta m`$. Srivastava et al.  derived a relation that differs by at least a factor of 2 from the standard result. The origin of this factor was studied by Lowe et al.  and by Burkhardt et al.  who found an error in Refs. , and demonstrated that the standard result is recovered when this error is corrected.

Other treatments  have studied some consequences of differing assumptions about the system and, in particular, the energies and momenta of the mass eigenstates. Since the $`K_L`$ and $`K_S`$ states that make up the oscillating $`K^0\overline{K}^0`$ system have different masses, they cannot have both the same momentum and the same energy. For example, Lipkin and collaborators  have studied the consequences of assuming either equal momentum or equal energy for the $`K_L`$ and $`K_S`$. In several of the above papers these two kinematic assumptions are examined, and some papers predict a wavelength for the oscillations which differs by a factor of exactly 2 from the standard treatment. A recent paper by Okun et al.  gives a concise summary of the situation.

However, we are not free to choose the energy and momentum of the mass eigenstates. Usually, the neutral particles are produced either in a reaction (e.g. $`\pi ^{}p\mathrm{\Lambda }K^0`$) for the case of kaons, or a decay (e.g. $`\pi \mu \nu _\mu `$) for neutrinos. In either case, the mass eigenstates have neither the same energy nor the same momentum for a given center-of-mass energy (mass in the rest frame of the source). The energies and momenta are determined by the kinematics. This was pointed out initially by Boehm and Vogel , Goldman , Srivastava et al. , Dolgov  and also by the present authors , who showed that when this is taken into account correctly, a consistent treatment of the kinematics follows and the standard result for the wavelength of oscillations is recovered.

However, Lipkin et al.  and Stodolsky  claim that the two mass eigenstates must have the same energy. In particular, Lipkin  states that mass eigenstates with different energy cannot interfere to produce the oscillations. In sect. 2 we show that this argument is incorrect and that the $`K_L`$ and $`K_S`$ states can indeed interfere to give strangeness oscillations. In sect. 3, we show that the mass eigenstates with kinematically correct energies and momenta produce oscillations with the same wavelength as in the standard treatment, without any additional factors.

## 2 Coherence of the mass eigenstates

In this section, we examine the coherence of two interfering wave functions. Suppose the wave functions are plane waves, $`\psi _1`$ and $`\psi _2`$. These might be, for example, two parts of the wave function in a 2-slit optical experiment, or an electron diffraction experiment. If so, they will have the same energy and momentum, but for generality, we keep both energies and momenta distinct for now.

The wave function at the point of interference, $`(x,t)`$, is

$`\psi =\psi _1+\psi _2=\sqrt{{\displaystyle \frac{1}{2}}}\left[\mathrm{exp}\{\mathrm{i}(p_1xE_1t)\}+\mathrm{exp}\{\mathrm{i}(p_2xE_2t+\varphi )\}\right]`$ (1)
where $`\varphi `$ is some phase angle introduced by the geometry (for example, the path difference between the two slits to the interference point or some difference induced at the presumed common source of the two components). We assume only that $`\varphi `$ is fixed and does not vary, e.g., randomly for different times due, for example, to fluctuations in the background medium (e.g., due to index of refraction fluctuations). It is in this sense that we may refer to the two waves as “coherent” . The probability density is

$`\psi ^2=1+\mathrm{cos}[(p_1p_2)x(E_1E_2)t+\varphi ]`$ (2)
In general, the second term oscillates in time with a time period characteristic of the energy scale of the system. For the optical case, this is $`10^{15}`$s, and in particle-physics experiments, the characteristic time is much shorter. In either case, this is well below the time resolution of any normal detector, so the second term averages to zero in the measurement, and $`\psi _1`$ and $`\psi _2`$ are therefore incoherent in the sense of Lipkin et al. . There may be other reasons for the cross term in Eq. (2) to vanish (e.g. different spin wave functions or different internal states of particles associated with $`\psi _1`$ and $`\psi _2`$), but in the absence of any such reason, the rapid time dependence of the energy term is the only reason why $`\psi _1`$ and $`\psi _2`$ are orthogonal when $`E_1`$ and $`E_2`$ are different. It is this energy term that seems to be the basis for statement often made (e.g. in ) that the kaon mass eigenstates are incoherent, and will not interfere, unless they have the same energy.

However, we must examine Eq. (2) in the particular case where $`\psi _1`$ and $`\psi _2`$ are $`K_L`$ and $`K_S`$ states. Here, $`E_1E_2`$ is of order $`\delta m`$, which is $`3\times 10^6`$eV. The associated time scale is $`2\times 10^{10}`$s, which is readily measurable. Thus, even in this plane-wave case (that is, in the absence of packeting), the response of a detector to $`\psi ^2`$ at a fixed point $`x`$ will oscillate measurably in time as given by Eq. (2). Furthermore, in actual experiments, the kaons do not appear in continuous plane waves, but as wave packets. Since kaon velocities in experiments are usually an appreciable fraction of $`c`$, the distance scale is many cm, again readily measurable. Of course, if the measurement averages over time or distance scales large compared with these values, then the cross term vanishes and the states become incoherent.

Thus, although the fact that the kinematics of the kaon (or neutrino) production process preclude equal energy for the two interfering states, the states may nevertheless be coherent, in the sense that we refer to above, and interference may be observable in certain situations.

## 3 Kinematics

Here, we calculate the kinematics for a specific case. For definiteness, we choose the $`K^0\overline{K}^0`$ example rather than neutrinos because

(i) There are just two states rather than three or more,

(ii) Accurate numerical values are known for the masses and the mass difference,

(iii) For neutrinos, it has been argued that a full treatment requires the inclusion of the detector in the system . However, for kaons, the $`\mathrm{\Delta }S=\mathrm{\Delta }Q`$ rule implies that the $`K^0`$ and $`\overline{K}^0`$ components can be identified from the kaon decay, without the need for a specific detector, thus simplifying the problem.

Neutrinos or kaons are produced either by a reaction such as

$$\pi p\mathrm{\Lambda }K^0$$

or by a decay, for example

$$\pi \mu \nu _\mu .$$

For the first of these, the center-of-mass energies and momenta of the mass eigenstates are given by

$`p_i^2={\displaystyle \frac{(sm_i^2m_\mathrm{\Lambda }^2)^24m_i^2m_\mathrm{\Lambda }^2}{4s}},E_i={\displaystyle \frac{s+m_i^2m_\mathrm{\Lambda }^2}{2\sqrt{s}}}`$ (3)
where $`i=S`$ or $`L`$. Similar expressions hold for the neutrinos from pion decay. The quantity $`\sqrt{s}`$ is the total center-of-mass energy for a reaction or the mass of the decaying particle for a decay. Although there may be a spread in the value of $`\sqrt{s}`$ for the overall system wave packet, our analysis proceeds component by component, i.e. at a precise value of $`\sqrt{s}`$ (within the constraints of the uncertainty principle). Thus $`p_Sp_L`$ and $`E_SE_L`$. In the following, to make the equations more readable, we ignore CP violation and we omit the widths of the kaon states.

Since the above reaction produces a pure $`K^0`$ state, the wave function at the reaction point, where $`x=t=0`$, is

$`K^0=\sqrt{{\displaystyle \frac{1}{2}}}(K_L+K_S).`$ (4)
This state develops in time as

$`\psi (x,t)=\sqrt{{\displaystyle \frac{1}{2}}}\{\mathrm{exp}[\mathrm{i}(p_LxE_Lt)]K_L+\mathrm{exp}[\mathrm{i}(p_SxE_St)]K_S\}.`$ (5)
Since Eqs. (3) give $`p_i`$ and $`E_i`$ in the center-of-mass frame, the $`x`$ and $`t`$ in Eq. (5) should also be in this frame. However, Eq. (5) and all following equations involve only invariants, so may be reinterpreted in the lab frame. Note that the two components in Eq. (5) are coherent in the sense that we defined above: no fluctuation of the relative phase occurs at the source. Under the assumption of propagation in a vacuum, there is also no medium to induce phase fluctuations coupled to the medium. At $`(x,t)`$, the probability amplitude for detecting a $`K^0`$ is

$`K^0\psi (x,t)={\displaystyle \frac{1}{2}}\{\mathrm{exp}[\mathrm{i}(p_LxE_Lt)]+\mathrm{exp}[\mathrm{i}(p_SxE_St)]\}`$ (6)
using $`K^0K_L=K^0K_S=\sqrt{\frac{1}{2}}`$. This is just as Eq. (1) above (with $`\varphi =0`$), giving the probability density as

$`K^0\psi (x,t)^2={\displaystyle \frac{1}{2}}\{1+\mathrm{cos}\left[(p_Lp_S)x(E_LE_S)t\right]\}.`$ (7)
Eq. (7) describes a plane-wave situation, with a unique value of $`\sqrt{s}`$, and with $`x`$ and $`t`$ as independent variables. In a realistic case, it would be used to form wave packets for all particles, with a spread of values of $`\sqrt{s}`$. The size of such wave packets might be determined, for example, by the time and position resolution of a detector in the incident beam or by the time structure of the accelerator beam. In any case, the packet cannot be larger than the $`K_S`$ lifetime, $`\tau _S0.9\times 10^{10}s`$, which gives a spacial extent of typically about 2 cm. Thus the outgoing kaon moves in a packet of this size, centered at the classical position. So the observation of the kaon at position $`x`$ must be made at a time when the packet is present, i.e. at time that is equal to $`x/\beta `$ or which differs from it by no more than the half-width of this wave packet. We therefore replace $`t`$ in Eq. (7) with the time defined in this way using beta calculated from the average $`K_L`$ and $`K_S`$ parameters:

$`\beta ={\displaystyle \frac{p_L+p_S}{E_L+E_S}},`$ (8)
which lies between the velocities $`p_L/E_L`$ and $`p_S/E_S`$. It is crucial that the observation is made at a single space-time point; no meaning can be attached to the interference of wave functions at different values of $`x`$ or $`t`$ (see ). In any realistic case, the separation of the centres of the $`K_L`$ and $`K_S`$ wave packets is much smaller than their size, so in practice, there is no loss of coherence due to separation of the $`K_L`$ and $`K_S`$ packets.

Thus Eq. (7) becomes

$`K^0\psi (x,t)^2={\displaystyle \frac{1}{2}}\left[1+\mathrm{cos}\{(p_Lp_S)x{\displaystyle \frac{(E_LE_S)}{\beta }}x\}\right].`$ (9)
Using

$`{\displaystyle \frac{1}{\beta }}={\displaystyle \frac{E_L^2E_S^2}{(p_L+p_S)(E_LE_S)}}={\displaystyle \frac{p_Lp_S}{E_LE_S}}+{\displaystyle \frac{2m\delta m}{(E_LE_S)(p_L+p_S)}},`$ (10)
Eq. (9) becomes

$`K^0\psi (x,t)^2={\displaystyle \frac{1}{2}}\left[1+\mathrm{cos}{\displaystyle \frac{m\delta m}{p}}x\right]`$ (11)
where $`m`$ and $`p`$ are the mean neutral kaon mass and momentum. Thus the wave number of strangeness oscillations in space is $`k=m\delta m/p`$. This is often expressed as an oscillation in time, in which case, the angular frequency is $`\omega =m\delta m/E`$. Both of these results are in agreement with the standard results , without any assumption about equality of momenta or energies.

## 4 Discussion

We conclude that a treatment of kaon and neutrino oscillations, with the kinematics treated in full, does indeed give the correct relation between the mass difference and the oscillation wavelength. Most of the recent literature is in agreement with this standard result for the wavelength, though sometimes using incorrect kinematics. If CP violation and the finite kaon lifetimes are incorporated in the above algebra, a more realistic result is obtained, with more cumbersome-looking equations, but the wavelength of the oscillations remains the same. The full equations are given in .

Lipkin  has suggested that the detector should also be included in the wave function of the system, since interaction with it depends on the details of the neutrino or kaon wave function. By choosing the kaon system here, we avoid this problem since the particles resulting from $`K^0`$ or $`\overline{K}^0`$ semileptonic decay identify the strangeness eigenstate, so no kaon detector, as such, is required. However, inclusion of a detector would not change our conclusions; at the zeros of the $`K^0`$ oscillation pattern, there are only $`\overline{K}^0`$ mesons and no $`K^0`$ mesons, so no detector could detect a $`K^0`$ at such a point.

If the mass eigenstates did in fact have equal momenta or equal energies, then this would imply a failure of 4-momentum conservation in the kaon or neutrino production process. Such a failure would be evident well outside the oscillation region; the $`K_L`$, which is the only state left in the asymptotic region, would, in principle, have an energy or momentum inconsistent with 4-momentum conservation.

In a recent preprint by Field , the kinematics and other aspects of the neutrino production are treated quite differently, giving a wide range of correction factors to the standard result for the wavelength. However, there is an error in his derivation of the pion decay rates (Eq. (7) of Ref. ). When this is corrected, the motivation for his later modifications to the standard treatment, with their bizarre physical consequences, is removed.

This research is supported in part by the Department of Energy under contracts W-7405-ENG-36 and DE-FG03-94ER40821, and by the NSF under contract PHY-0099385.


text:
# Interplay between geometry and flow distribution in an airway tree

## Abstract

Uniform flow distribution in a symmetric volume can be realized through a symmetric branched tree. It is shown here however, by 3D numerical simulation of the Navier-Stokes equations, that the flow partitioning can be highly sensitive to deviations from exact symmetry if inertial effects are present. The flow asymmetry is quantified and found to depend on the Reynolds number. Moreover, for a given Reynolds number, we show that the flow distribution depends on the aspect ratio of the branching elements as well as their angular arrangement. Our results indicate that physiological variability should be severely restricted in order to ensure adequate fluid distribution through a tree.

The problem of fluid flow in a branching geometry appears in many physical, geological, chemical and biological systems. Examples include catalysis, flow through porous media, blood circulation and respiration. When studying transport in any of these systems, a common objective is to understand the mechanisms that govern the flow partitioning at the interconnections level. Until recently, it has been generally assumed that the use of Darcy’s law should be sufficient to describe the propagation of flow through branched structures. Such a relation corresponds to the linear dependence between flow and pressure drop, $`Q\mathrm{\Delta }P`$, which is strictly valid at small Reynolds number. Regardless of this limitation, a large number of studies have been based on this approximation. In the context of porous media, for instance, a simple paradigm to represent flow through the pore space is a network of bifurcating and merging channels where the transport of fluid is analogous to the distribution of electrical currents in a resistor network. However these models can only predict a perfectly uniform and synchronous flow distribution through airways bifurcations . A major problem in modeling of flow through trees arises from the fact that, due to inertial effects, Darcy’s law breaks down as a phenomenological description for large Reynolds numbers. Even at moderate Reynolds, the inertial non-linearities become relevant as compared to viscous effects.

Unambiguous experimental and numerical evidences of inertial effects have been observed in several studies on flow though branched structures, with special emphasis on the bronchial tree . Such phenomena exists in real lungs but they are more simple to study in a symmetric geometry . In particular, in order to irrigate uniformly a symmetric volume it is easy to show, through the following collage argument, that this is ensured by a symmetric tree. Suppose that an asymmetric tree feeds a volume which has a plane of symmetry. If the tree is asymmetric the flow will be different in the two parts of the volume which are symmetrical. Then one can replace the tree with a non-uniform flow by the symmetry image of the more efficient region. The new tree, which is now symmetric, is more efficient for flow distribution.

In the Poiseuille approximation, the only way to have perfect symmetry is to work with an equivalent resistor network that is symmetric. In other words, at each bifurcation the daughter branches should be exactly identical irrespectively of their real geometrical arrangement. This might not be true if inertial effects are present. It should be recalled that, as the lung is a succession of branch bifurcations, the final flow distribution can be represented by a multiplicative process. In consequence, even a rather small asymmetry could lead to a strong inhomogeneity of the flow distribution . Because the geometrical arrangement of the bronchial tree of mammals is always subjected to some physiological variability , it appears natural to question whether a small modification of the structure disturbs the distribution of fluid flow.

The purpose of this work is to investigate how the tree geometry influences the flow in order to shed some light on the optimal aspect of the bronchial tree for distributing air uniformly in the lung volume. The direct 3D numerical solution of the Navier-Stokes equations is by far the most practical way to elucidate this problem. The simplified tree model used here is shown in Fig. 1. It consists in a 3-dimensional cascade of cylinders branching through two bifurcations. Each bifurcation ABC or BDE or CFG, is coplanar as found approximately in real lungs. The bifurcation geometries are modeled in such a way to minimize geometrical singularities as shown in Fig. 2. For simplicity, we assume that the radii of the tubes decrease with a factor $`2^{1/3}`$ at each bifurcation  and choose the branching angle to be $`45^{}`$.

The mathematical description for the detailed fluid mechanics in the branched structure is based on the steady-state form of the continuity and Navier-Stokes equations for mass and momentum conservation

$$𝐮=0,$$
(1)
$$\rho 𝐮𝐮=p+\mu ^2𝐮,$$
(2)
where u and $`p`$ are the local velocity and pressure fields, respectively. The diameter of the first tube is equal to $`2`$ cm, corresponding approximately to the diameter of the human trachea. The fluid is air with viscosity $`\mu =1.785\times 10^5`$ $`\mathrm{kg}\mathrm{m}^1\mathrm{s}^1`$ and density $`\rho =1.18`$ $`\mathrm{kg}\mathrm{m}^3`$, and the flow is considered to be incompressible. Nonslip boundary conditions are imposed at the tube walls (Dirichlet condition $`𝐮=0`$) and the velocity at the entrance A is parabolic. The outlets are free with the same reference pressure and $`𝐮/n=0`$. Equations (1) and (2) are solved using finite elements . For all simulations, the relative conservation error is smaller than $`3\%`$.

The parameters governing the flow are the bronchi aspect ratio (length to diameter ratio of the tubes) $`L/D`$, the rotation angle $`\alpha `$ between successive bifurcations, and the Reynolds number, $`Re\rho DV/\mu `$, where $`V`$ is the mean velocity at the entrance. The reference angle $`\alpha =0^{}`$ corresponds to a coplanar tree. The flow asymmetry is defined as

$$\mathrm{\Sigma }(\alpha ,L/D)\left|\frac{q_1q_2}{q_1+q_2}\right|,$$
(3)
where $`q_1`$ and $`q_2`$ are the outflows at (D,G) and (E,F) branches, respectively. We perform simulations for several values of $`\alpha `$, $`L/D`$ and $`Re`$ to find their influence on the flow partitioning $`\mathrm{\Sigma }`$. Note that the air velocity at the entrance of human lungs range from $`1`$ m/s at rest ($`Re1200`$) to $`10`$ m/s for the condition of very hard exercise ($`Re12000`$) . Due to the number of parameters governing the flow and the computation time for each set of parameters, we first discuss the dependence of the flow asymmetry on the geometry for a fixed Reynolds value, namely, $`Re=1200`$. This corresponds to the human inspiration state at rest.

The results are shown in Fig. 3. The main result is that, whatever the conditions, the behavior of $`\mathrm{\Sigma }`$ around the minimum is not parabolic. Even a small departure from geometrical symmetry can cause a non negligible flow disturbance. For a given value of $`L/D`$, the disturbance increases with the deviation of the rotation angle from $`90^{}`$. $`\mathrm{\Sigma }`$ is therefore maximum for a planar tree and, for a fixed $`\alpha `$ value, it decreases with increasing aspect ratio.

There are then two facts to interpret. First, why the flow is influenced by breaking the symmetry only. Seconds why this effect is attenuated for long branches or large aspect ratios. The first fact can be understood by considering the velocity distribution in a cut of the secondary branch B shown in Fig. 4. The flow keeps the symmetry of the ABC bifurcation plane but, due to inertia, the high velocity regions are drifted vertically and an M-shape type of distribution is observed . This shape governs the flow partitioning at the second bifurcation. Note that if the branches B and C are long enough and for small $`Re`$, the profile should tend to a parabolic type. As a consequence, the distribution shown in Fig. 3 will progressively change along the second generation branch. It is because the branch length is too short that the grand-daughter branches can capture the asymmetry seen in Fig. 4. This provides a qualitative answer for the second question. The position of the intersection relatively to the M-shape is then the key for asymmetry. For example in Fig. 4, the branch E obviously receives more flow than branch D. It is also clear that, if $`\alpha =90^{}`$, the flow symmetry is restored for any value of $`L/D`$.

The dependence of the flow asymmetry as a function of the Reynolds number is shown in Fig. 5. A strong increase of $`\mathrm{\Sigma }`$ is observed up to $`Re250`$ followed by a region of weaker dependence. This type of behavior has been previously reported for 2D flow in trees comprising more than two generations of branches . It is remarkable that the onset at $`Re250`$ is approximately the same whatever the angle $`\alpha `$. This is a clear indication that, at the entrance of the second bifurcation, the velocity profile reaches the same pattern for a given $`Re`$ value. Again, the final asymmetry of the distribution of flow is a result of the inertial effects originated from the first bifurcation. All these arguments are illustrated in Fig. 6, where the contour plots of the velocity fields are shown at the entrance of the second bifurcation. At large $`Re`$, the M-shape is revealed and, as expected, the lower the $`Re`$, the closer the profile is to parabolic flow. The smaller variation of $`\mathrm{\Sigma }`$ for $`Re>250`$ can be explained by the presence of a secondary flow .

Some implications of our results are noteworthy. If the inertial effects observed here are present in a larger tree, the relative flows delivered to the outlets of this structure may become strongly non-uniform. This broadness in the flow distribution is a typical signature of a multiplicative process , where an observable can be viewed as a “grand process” depending on the successful completion of a number $`n`$ of independent “subprocesses”. It is then possible to associate the flow at each branch with a probability $`p_i`$, so that the flow at a given outlet $`k`$ be $`q_kp_1p_2\mathrm{}p_n`$, where $`i=1,2,3,\mathrm{},n`$ corresponds to the set of branches constituting the pathway going from the entrance to the exit $`k`$. It can be easily demonstrated that, if the $`p_i`$’s are independent variables and $`n`$ large, the distribution of $`q_k`$ should be approximately log-normal. Furthermore, this distribution might mimic a power-law if its dispersion is sufficiently large . Note that this situation is that of the human bronchial tree (where $`L/D`$ is close to $`3`$) even at rest. In this case $`Re1200`$ and the multiplicative process due to inertia can propagate further down in the tree. If we consider that these effects only disappear for $`Re`$ less than 100 and that the local Reynolds decreases by a factor of $`2^{2/3}`$ at each generation, we obtain that the flow asymmetry can be significant up to the 6th generation of the bronchial tree under rest conditions .

In conclusion, we have investigated the effect of inertia on fluid flow through three-dimensional rigid branched structures by direct numerical simulation of the Navier-Stokes equations. It has been found that for trees with $`3`$ generations of cylindrical conduits, the flow distribution at the outlets strongly depends on the Reynolds number and on the geometry of the ramified structure. Moreover, our simulations indicate that the flow imbalance throughout the tree is highly sensitive to the aspect ratio $`L/D`$ of its cylindrical units and to the variation of the rotation angle $`\alpha `$ between successive bifurcations. While a uniform distribution of flows at the outlets of the third generation branches is always obtained for $`\alpha =90^{}`$, our calculations show that a small deviation from this geometrical configuration is capable to induce a large asymmetry on the flow. Note that the presence of long branches would lead to purely axi-symmetric parabolic profiles and flow symmetry. However, long tubes exhibit large hydrodynamic resistance (proportional to $`L/D^4`$). It is therefore not surprising that in real lungs $`L/D3`$ and $`\alpha 90^{}`$.

Finally, our results suggest that small deviations from the “best” structure should have the same type of consequences in the real (asymmetric) lung, namely, strong dependence on geometry and Reynolds number. In particular, the flow distribution at rest and exercise might be significantly different. These results could also help to better understand lung morphology. It has been argued  that the asymmetric structure of the lung is solely due to geometrical constraints. Our study indicates that the inertial effects plays also an important role in air distribution. In other words, the asymmetry of the bronchial tree is determined not only by geometrical constraints but also by the existence of inertial effects. Of course, if the flow distribution is found uniform although the geometry is “imperfect”, the following question would naturally arise : what are the physiological regulation mechanisms that can compensate the flow non-uniformity due to inertial effects? In addition, the fluid dynamics studied here is certainly relevant to understand particle deposition in the airway tree , a problem of crucial importance both from the physiologic and the therapeutic points of view.

We thank CNPq, CAPES, COFECUB and FUNCAP for support. The Centre de Mathématiques et de leurs Applications and the Laboratoire de Physique de la Matière Condensée are “Unité Mixte de Recherches du Centre National de la Recherche Scientifique” no. 8536 and 7643.