Tanyu Bonev
email: tbonev@astro.bas.bg

Klaus Jockers
email: jockers@mps.mpg.de

Nikolay Karpov
email: karpov@mao.kiev.ua

Klaus Bonev

Nikolay Jockers

Karpov

A dynamical model with a new inversion technique applied to observations of comet C/2000 WM 1 (LINEAR)

Keywords: COMETS, COMETS, DUST, COMETS, DYNAMICS

Three continuum images of comet C/2000 WM 1 (LINEAR) obtained on Nov 10, Nov 19, and Dec 03, 2001, are analyzed with the aid of a dynamical model, i.e. with a model that uses the size-dependent motion of dust grains under solar radiation pressure to determine the dust size distribution and its temporal change. The frames are photometrically calibrated in terms of the albedo

filling factor product. On Nov 20.2 the Earth transited the orbital plane of the comet and an anti-tail was recognized in the image of Nov 19. For the determination of the particle fluxes describing the contribution of monodisperse particle shells to the cometary brightness the model uses a new regularization method employing Chebyshev polynomials of selected orders in emission time and particle size. It guarantees positiveness of the particle fluxes and imposes a varying degree of smoothness on their dependence on particle size and emission time. The particle emission velocities are still derived by trial and error. The dynamical model is described in detail. Results are presented for several low orders of the Chebyshev polynomials and are compared in order to understand the limitations imposed by the regularization process. The size distributions derived from the different observations do not always agree. This is particularly true for the earliest and most recent synchrones contributing to an image. In the observations of Nov 10 and Dec 03, i.e. excluding the anti-tail image, the integrated mass loss strongly decreases in the most recent time steps of the model although the comet is still approaching the Sun. This is interpreted as an artifact introduced by the overlap of the shells of large particle size emitted shortly before the observation. The model derives an increasing number of small particles released by the comet in the second half of November. This is at least in part considered as real and attributed to parti-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 1 Introduction

Dust particles in the coma and tail of comets scatter the solar light and give rise to the cometary continuum at optical wavelengths. Induced by the scattering process the particles experience a radiation pressure force depending on their linear dimension, repelling them from the Sun. Dynamical models of the surface brightness distribution of cometary dust make use of this fact and, if the dust particle density and albedo are known, can provide information on the number and mass of dust particles released from the nucleus and their size distribution as a function of time. For the dynamical modeling it is usually assumed that the dust particles remain unchanged during their motion in the cometary coma, i. e. do not release volatile material and do not fractionate.

The work by Finson and Probstein (1968b); [START_REF] Kimura | On the structure of cometary dust tails. Chinese A C C E P T E D M A N U S C R I[END_REF]; [START_REF] Liu | The Structure of Dust Tails of Comets -Part Two -the Tail and Dust Content of Comet Arend-Roland[END_REF] and [START_REF] Jewitt | A CCD portrait of Comet P/Tempel 2[END_REF] provides early examples of successful application of dynamical modeling for the description of cometary dust comae and tails.

In the optically thin dust coma particles of different sizes ejected at different times contribute linearly to an observed image:

b = M • F (1)
b is the observed brightness of a cometary dust image and M a vector of model images. The model images are images of particle clouds ("shells") valid at the time of the observed cometary dust image. Each shell belongs to a certain emission time and particle size. In most models their ejection velocity is uniquely determined by their particle size. F is the flux vector that determines how many particles of a certain size and ejection time are present in the 1) nonunique. In the approach of direct modeling plausible assumptions are made about the grain size distribution and its temporal evolution. The derived models are then compared with the observed brightness distribution and, if necessary, corrected, until satisfactory agreement between model and observation is achieved.

A major step forward in the analysis of images of cometary dust comae and tails occurred when Fulle (1987Fulle ( , 1989) ) applied the method of [START_REF] Tikhonov | Sublimating components in the coma of comet C/2000 WM 1 (LINEAR)[END_REF] and solved the inverse problem by regularization. The regularization technique is supposed to find the "most probable" flux vector, i.

e. the "most probable" particle size distribution and its temporal evolution.

In any case, because of its non-linear relation to the flux vector, the particle ejection speed, its dependence on particle mass and a possible anisotropy of the ejection of particles from the nucleus must still be derived by trial and error. Fulle has applied his method successfully to a large number of comets (for an overview see Fulle, 2004).

Despite the great success of Fulle's inverse dynamical modeling no publication of other authors using his regularization method has appeared in the literature 1) by trial and error or by other methods (for more details see Section 5.5). In this paper, we employ for the inversion of eq. ( 1) a two-dimensional generalization (for the two quantities emission time and particle size) of the one-dimensional technique developed by [START_REF] Lemen | Coronal activity in F-, G-, and K-type stars. III -The coronal differential emission measure distribution of Capella, Sigma-squared CrB, and Procyon[END_REF] for derivation of plasma parameters from comparison of observed and theoretical X-ray spectra of the coronae of cool stars. Like in other models applied previously, in our model the particle ejection speed and its dependence on particle mass is derived by trial and error. Lemen's method allows to enforce a varying degree of smoothness on the dependence of the flux vector F on time and particle size by employing Chebyshev polynomials of different order for the fits. In this way the method allows to derive several solutions of the inverse problem based on different orders in time and particle size. Comparison of such solutions provides insight into the accuracy of the inverse solutions. Therefore, for the first time we can investigate "what the regularization does to the model" (Jockers, 1997b).

The long-period comet C/2000 WM 1 (LINEAR) (in the following comet WM 1 ) was discovered by the Lincoln Near Earth Asteroid Research (LINEAR) facility on Dec 16, 2000. The comet reached its perihelion of 0.555 AU on Jan 22.6734, 2002 [START_REF] Marsden | The Temporal Development of the Pre-perihelion Infrared Spectral A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT Energy Distribution of Comet Hyakutake (C/1996 B2)[END_REF]. The comet was favorably placed in the sky for northern hemisphere observers in November and early December 2001 when it was still outside of the Earth's orbit approaching perihelion. It reached its minimum geocentric distance of 0.316 AU on December 2.6, 2001.

On Nov 20.20 UT the Earth passed through the orbit plane of the comet. A broad anti-tail consisting of old heavy particles, trailing behind the comet, was visible. [START_REF] Watanabe | Comet C/2000 WM 1 (LINEAR)[END_REF] has reported lengths and position angles of this anti-tail. His results are based on I-band CCD images obtained in the period In Section 2 we briefly describe the observations and present the three input images calibrated in dimensionless units of the albedo filling factor product (Figure 1). Then (Section 3) we calculate the synchrones and syndynes belonging to the three observations and use them for a first image analysis (Figure 3). In Section 4 we describe our model. 4).

In Section 5.6 the appearance of a large number of small particles in comet WM 1 in the second half of November is put in context with other observations of WM 1 and with similar events in other comets (Figure 15). Section 6 summarizes and concludes the paper.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

The Observations

The observations presented in this paper were carried out with the 2-m telescope of Pik Terskol Observatory. The Two-Channel Focal Reducer (Jockers, 1997a) a size of 0.027 mm and therefore provides a scale of 0.99 arcsec px -1 .

The three observations were obtained 73, 65, and 50 days before perihelion, respectively. Table 1 gives a summary of the observing circumstances.

The images were flat fielded and the sky brightness, measured at image areas far from the comet, was subtracted. Star trails were marked and their areas were replaced by adjacent areas of similar count rate, but free of stars. The photometric calibration was performed with the help of calibration images of spectrophotometric standard stars from the list of [START_REF] Hamuy | Southern spectrophotometric standards[END_REF][START_REF] Hamuy | Southern spectrophotometric standards, 2[END_REF] and [START_REF] Bessell | Spectrophotometry: Revised Standards and Techniques[END_REF] taken at different airmasses during the observations. In a next step the images were expressed in units of the albedo filling factor product pΦ(φ)f (x, y). In this product p denotes the geometric albedo of the dust particle, i.e. the albedo at phase angle φ = 0, Φ(φ) is the phase function at phase angle φ, normalized to unity at φ = 0. f (x, y) denotes the local filling factor, i.e. the fraction of sky area covered by the geometric cross section of the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
dust particles. The albedo filling factor product pΦ(φ)f (x, y) is connected with the intensity i(x, y) of the calibrated images through the equation (Jockers, 1997b):

p Φ(φ) f (x, y) = i(x, y) i r h R 2 . ( 2
)
r h is the heliocentric distance, R and i are the radius and the mean disk intensity of the Sun, respectively [START_REF] Allen | Astrophysical quantities[END_REF]. Figure 1

Syndynes-synchrones analysis

We start the dynamical analysis with the zero ejection velocity approximation and calculate the syndynes and synchrones of the dust particles using the equations given by Finson and Probstein (1968a). After being released from the nucleus, particles move in the gravitational field of the Sun (gravitational force F g ) reduced by the force of radiation pressure F r . The effective force

F A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
exerted on a dust particle is given by

F = F g -F r = F g • (1 -β), (3)
where the ratio of radiation pressure to gravitational force, β, can be calculated from [START_REF] Burns | Radiation forces on small particles in the solar system[END_REF]:

β = 1.19 10 -4 • Q pr /(d • ρ). ( 4
)
Here Q pr is the radiation pressure efficiency, and d is the diameter of the sphere, expressed in cm, having the same volume as the dust particle. ρ is the particle density in g cm -3 . In this paper we adopt ρ = 1 g cm -3 . Provided a β value has been derived from observations and Q pr is known from theoretical considerations, the product of the diameter and density of the particle can be determined. An empirical radiation pressure efficiency, valid in a wide size range, was derived by Bonev and Jockers (2002). These authors used oblate spheroids with axial ratio 1.8 as test particles. The optical properties of the particles were characterized by a refractive index equal to 1.52 + i × 0.01. For the small particles (0.2μm < d < 3.2μm), the calculation of Q pr was performed with the T-matrix method. Geometrical optics calculus was applied in the case of large particles (12μm < d < 90μm). As it is still impossible to calculate the optical properties of particles of intermediate size, Q pr was interpolated in this size range. The results of these calculations are well reproduced by the following equation (Bonev and Jockers, 2002):

Q pr = a 3 1 - a 0 d a 1 a 0 d a 2 p(a 1 , a 2 ) + a 4 d a 5 + d p(a 1 , a 2 ) = 1 + a 1 a 2 a 2 1 + a 2 a 1 a 1 , d ≥ a 0 . (5) A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
In eq. ( 5) the parameters have the following values: a 0 = 0.02 μm, a 1 = 34.848, a 2 = 1.094, a 3 = 0.61, a 4 = 0.7, and a 5 = 2.0 μm. Note that some of the parameters differ from the value provided by Bonev et al. (2002), as we deal here with diameters, not with radii. is the one starting from the cometary nucleus closest to the anti-solar direction.

The corresponding particle diameters (eqs. 4, 5) are d = 5.0, 12.5, 32, 82, 207, and 520 μm. The radiation pressure efficiency Q pr varies between 0.699 and 0.669, i.e. very little. The synchrones (continuous lines) delimit the synchrone ranges we will use in section 4 in the Monte Carlo modeling. The numbers between the synchrones refer to Table 2 where they are listed together with the start and end times of the synchrone ranges and the corresponding true anomaly and heliocentric distance ranges of the comet.

On Nov 10 the shape of the tail is well described by the syndynes with βvalues 0.025, 0.010, and 0.004, (particle sizes 32 -207 μm). Further down, closer to the anti-solar direction (the Sun is to the left), smaller (μm-sized) and younger particles, released from the nucleus a fraction of one day to several days before the observation, are expected. From the decrease of the pΦf values in this region we conclude that on Nov 10 the comet should be characterized by a deficiency of small particles.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

On Nov 19 the most remarkable feature is the anti-tail. It is extended close to the projected direction to the Sun and along the negative velocity vector of the comet. The shape of the anti-tail is well covered by the same three syndynes as on Nov 10. Possibly even larger particles with smaller β-values and syndynes closer to the projection of the cometary orbit are needed to accurately describe the details in this image. On Nov 19 the Earth was only 2.5 degree below the orbital plane of the comet. The syndyne analysis shows that the observed anti-tail comprises large particles, located outside of its orbit close to the cometary orbit plane and trailing behind the cometary nucleus with respect to its orbital motion. Thus, the observed anti-tail can possibly be considered as the inner part of a trail, a term usually used to describe the distribution of very large particles trailing behind the nucleus of short period comets [START_REF] Sykes | The interplanetary dust complex and comets[END_REF][START_REF] Reach | Bayesian-based iterative method of image restoration[END_REF][START_REF] Agarwal | Dust Environment Modelling of Comet 67P/Churyumov-Gerasimenko[END_REF]. Here we observe it in a long period comet. According to Bredichin's terminology the observed elongation is a "pseudo" anomalous tail (see Festou et al. (1993); [START_REF] Jaegermann | Bredichins mechanische Untersuchungen üeber Kometenformen in systematischer Darstellung[END_REF]). The orientation of the anti-tail, derived from our image, coincides with the orientation of the anti-tail in the I-band images obtained by [START_REF] Watanabe | Comet C/2000 WM 1 (LINEAR)[END_REF]. From Table 2 we see that the oldest and heaviest particles have been released from the nucleus more than 290 days before perihelion, i.e. more than 230 days before the observation. Particles as large as 200 μm must be included in the Monte Carlo models.

In contrast to the images from Nov 10 and Nov 19 the image of Dec 03 exhibits a well expressed anti-solar, normal dust tail. This tail coincides well with the syndyne with β-value 0.16 (particle size 5 μm, significantly smaller than in the previous images). The release age of these particles does not exceed 15 days (Table 2). Such small, high β particles would populate the anti-solar direction.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

But the intensity distribution on Dec 03 shows an anti-solar elongation only close to the nucleus. This reflects merely the circumstance that all syndynes leave the nucleus exactly in anti-solar direction (Finson and Probstein, 1968a).

Further tailward they are bent away from this direction. This shape of the tail suggests again that the number of small young particles (diameter The question if the amount of large particles is indeed reduced in the Dec 3 observation as compared to Nov 19 must be addressed by more elaborate analysis, as will be done with our inverse Monte Carlo model.

<

Monte Carlo modeling

Direct Monte Carlo dynamical simulation

The direct Monte Carlo model calculates a synthetic image of the dust coma and/or tail of a comet for comparison with telescopic observed images. Its two major elements are the modeling of the ejection of dust grains from the cometary nucleus and the calculation of their motion under gravity reduced by solar radiation pressure.

A C C E P T E D M

A N U S C R I P T ACCEPTED MANUSCRIPT

Ejection of particles from the nucleus

The ejection of particles from the nucleus, described mathematically as the specification of the velocity vector with which the particles are released from the nucleus, is a rather complicated process, depending on many parameters.

On one hand we want to use the model to extract information about a real comet. Therefore it seems that one should describe the particle ejection from the nucleus as detailed as possible. On the other hand it is well known that the inverse problem is improperly posed. Therefore we should restrict the input parameters as much as possible to avoid that our model will become just a multi-parameter fit and that the derived parameters loose their physical meaning. In our case the model allows isotropic emission and emission into a cone with a prescribed opening angle and axis. In the present work one half of the particles were emitted isotropically into all directions, and the other half evenly into a cone of 30 • half opening angle directed to the Sun.

For both, the isotropic emission as well as the emission into a cone the ejection speed is described by the following equation:

v(d, r h , θ) = u(r h ) f (d) h(θ) ( 6 )
r h is the heliocentric distance, d the volume equivalent diameter of the particle, and θ the solar zenith angle at the place of ejection on the nucleus. [START_REF] Waniak | Nuclear dust emission pattern of Comet Wilson 1987VII[END_REF] has developed an algorithm based on the maximum likelihood approach that allows to derive the directional distribution of the dust emission and applied it to images of comet C/1986 P1 (Wilson). As there was no obvious indication for asymmetric dust emission in our images of WM 1 we assume the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
dependence of the emission velocity on the angle θ as simple as follows:

h(θ) = cos(θ/2). (7)
This functional form allows a finite emission velocity on the night side and has a squareroot singularity of zero ejection speed in anti-solar direction. In principle a zero outflow speed at a finite dust production rate produces a singularity in the dust pattern close to the nucleus. But, as the singularity in our case is of squareroot character and occurs only in the part of the distribution emitted isotropically in all directions, the effect is weak. It is further reduced because the radiative acceleration quickly removes even particles of zero emission speed from the neighborhood of the cometary nucleus.

The function f (d) describes how the ejection speed depends on particle diameter. Finson and Probstein (1968a) realized that the dust particles are accelerated by gas drag and therefore their final velocities depend on the gas flux at the nucleus. As in turn the gas flux is affected by the mass loading caused by the presence of the dust grains, the dust ejection speed depends on the quantity and distribution of the dust grains as well. In the framework of dynamical modelling the grain ejection speed must be uniquely determined by its size and we must neglect the feedback of the grains on the gas flow.

Therefore our approach is strictly valid only for a comet with low dust/gas ratio. In the framework of our model we describe the dependence of the ejection velocity on particle size f (d) by the following equation:

f (d) = 1.0/(1. + (d/d 0 ) κ ). ( 8
)
The parameter d 0 is the transition diameter from constant velocities for the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
smallest particles to a power law dependence for the larger ones, i.e. parameterizes the gas flow. In this paper d 0 = 0.25 μm. We found that for κ = 0.5 eq. ( 8) well fits the velocities derived by [START_REF] Crifo | Hydrodynamic models of the collisional coma[END_REF] and [START_REF] Gombosi | A heuristic model of the Comet Halley dust size distribution[END_REF] in their models of the gas-dust interaction in the circumnuclear coma. At the same time the relation is consistent with the approximation of [START_REF] Liu | The Structure of Dust Tails of Comets -Part Two -the Tail and Dust Content of Comet Arend-Roland[END_REF] and allows to treat cases of different dust-to-gas ratios and different gas production rates (Bonev and Jockers, 2002) relatively simply. Figure 5 shows the dependence given by eq. ( 8) for two different values of the dust-to-gas ratio. Figure 5 contains also the radiation pressure acceleration β.

For the specific case of our observations the dependence of the ejection speed on heliocentric distance u(r h ) was determined by trial and error. Trial shells were calculated for different particle sizes and emission times (which in our preperihelion observations are uniquely related to heliocentric distance). Both dependencies are well decoupled: The location of the centers of the shells is given by the synchrones corresponding to their emission time and heliocentric distance, and their position on the synchrone is related to their particle size.

The empirical velocity was selected to minimize overlap of the resulting shells and to provide optimum coverage of the observed dust coma as a whole. The empirical velocities are plotted in Fig. 4 (full line). They can be approximated by a power law with exponent -1.5:

u(r h ) = v max (r h /1.0 AU) -1.5 (9)
In this equation the maximum velocity v max , i. e. the maximum velocity at 1 AU heliocentric distance, = 500 m s -1 . The r -1.5 h dependence lies between the r -1 h dependence suggested by [START_REF] Whipple | A Comet Model. II. Physical Relations for Comets and Meteors[END_REF] and the r -2 h derived by Crifo and Rodionov (1997) (see discussion on the subject in Fulle ( 2000)). Our

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
data, in particular the anti-tail image, contain particles emitted at heliocentric distances > 4 AU. These particles restrict the exponent in heliocentric distance particularly strongly. The power law velocity is shown as dashed line in Fig. 2 below shows in columns 8 the quantity (r h /1.0 AU) -1.5 , in column 9 (v 1 ) the velocity according to eq. ( 9) and in column 10 (v 2 ) the empirical values found by trial and error.

Table

As the initial velocity vector enters the model in a non-linear way it must always be determined by trial and error. Inverse models (including ours)

presently cannot produce the initial velocity vector and its dependence on particle radius, gas to dust ratio, and heliocentric distance.

Particle shells and quasisynchrones

To calculate a model image we must emit particle clouds from the cometary nucleus at a number of time steps earlier than the time "of observation" of the model image t obs and calculate where these particles are at that time. The emission times are in principle arbitrary, one could, e. g., select them randomly. In our model we select the emission times such that they correspond to equal intervals of true anomaly. According to the second law of Kepler, if a constant number of particles is emitted at equidistant steps of true anomaly the production rate will vary as the inverse square of the heliocentric distance r -2 h . Such a variation is expected if all solar radiation received by the nucleus is expended for sublimation of water as should be the case at heliocentric distances < 2 AU, when the temperature of radiative equilibrium well exceeds the sublimation temperature of water. Emitting at time intervals equally spaced in true anomaly does not restrict modeling of outbursts, as the number of

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
emitted particles, as determined later by the inverse model, can vary at each time step.

In an approach similar to that of [START_REF] Moreno | Dust in Comet 67P/Churyumov-Gerasimenko[END_REF] we divide the emission times of particles present in an observed image into discrete time intervals.

We will call the dust particles emitted in a time interval from t q to t q + Δt q and "observed" at the time t obs ≥ t q + Δt q a quasisynchrone. In contrast to a classical synchrone, in a quasisynchrone the particles are emitted in a finite time interval instead of a single time instant and the initial speed is not zero but determined by eqs. ( 6) to (7) (instead of eq. ( 9) the empirical velocities of Fig. 4 are used). In our paper a quasisynchrone contains 100 time steps. As opposed to [START_REF] Moreno | Dust in Comet 67P/Churyumov-Gerasimenko[END_REF], in order to implement the r Figure 3). The β values are the same for each of the 100 time steps selected within the interval, actually they are the same for all quasisynchrones. In a quasisynchrone the particle size distribution (i. e. the numbers with which the shells belonging to a particle size bin are multiplied) is specified according to a power law. It is the same for all time steps within the quasisynchrone. Later the factors with which the shells must be multiplied will be the result of the inverse calculation.

In our models of comet WM 1 11 quasisynchrones are emitted in total. An overview on these quasisynchrones is provided in Table 2. The true anomaly

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
ranges (columns 2 and 3) are selected equidistant (excepting the intervals of the three observed images). The corresponding times and heliocentric distances are provided in columns 4 and 5, and 6 and 7, respectively. Column 8

indicates the values of r -1.5 H /1.0 AU. The velocity v 1 = v max r -1.5 H describes the smooth heliocentric dependence of outflow velocity and the velocity v 2 in column 9 the empirically found relation (see Fig. 4 and the comments following eq. ( 9)). For the three observed images n q = 8, 9, and 11 quasisynchrones are used, respectively.

In order to model one elementary shell i of a quasisynchrone q with particle size d i , at each of the 100 time subintervals 200 particles were emitted. Therefore within a quasisynchrone each shell with particles of diameter d i consists of The coordinates x and y have the dimension of length. n(x, y, t q , d i , t obs ) denotes the number of particles of the ith shell of the qth quasisynchrone found at observation time t obs in the projected detector bin with coordinates x and y.

N m = 100 ×200 = 2 •
We now can write down the model brightness of a particle shell with particle diameter d i emitted by quasisynchrone q and observed at time t obs :

M(x, y, t q , d i , t obs ) = C norm pΦ(φ) πd 2 i 4 g(t q , d i ) Δd i Δx Δy n(x, y, t q , d i , t obs ) (10)
pΦ(φ) is the product of geometric albedo and phase function, to be specified in sections 4.1.3 and 4.3.2 (see also Table 4). g(t, d) denotes the parti-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
cle size distribution function. It depends on time and is normalized so that g(t, d) dd = 1 for all t when integrated over all N s size shells with respect to their corresponding widths. For convenient comparison with observations it is useful to choose the normalization constant C norm = Δt q /N m , the inverse of the model production rate of particles of diameter d i . Then we get

M(x, y; t q , d i , t obs ) = pΦ(φ) πd 2 i 4 g(t q , d i ) Δd i Δx Δy n(x, y, t q , d i , t obs ) Δt q N m (11)
In the following we will refer to the quantity M(x, y; t q , d i , t obs ) as the normalized model brightness of a shell.

In total the final model surface brightness is formed as linear superposition of N q (8-11) quasisynchrones which in turn consist of a linear combination of 51 shells. Besides of the specification of the initial velocity discussed in section 4.1.1, which influences the final surface brightness numbers in a non-linear way and always must be found by trial and error, each quasisynchrone has a specific size distribution g(t q , d i ) consisting of N s numbers each. In the direct model the size distributions were specified as power laws g(t q , d i ) ∝ d ν . By

Albedo and phase function

The geometric albedo p is unknown. It cannot be determined from dynamical modeling as the cometary brightness depends on the product of particle number (filling factor) and albedo (eq. ( 2)). We adopted a value of 0.04 derived for the nucleus of comet Halley. No dependence of the albedo on particle size is assumed. Also the phase function is considered not to depend on particle size. This is justified by our preliminary analysis which has shown that only relatively large particles are contributing to the observed surface brightness of comet WM 1 . In accordance with [START_REF] Meech | Observations of Comet P/Halley at minimum phase angle[END_REF] we use a phase function of the form Φ(φ) = 10 -0.4 * γ * φ , where γ is the linear phase coefficient.

In the direct models we used γ = 0.035, the value derived by [START_REF] Meech | Observations of Comet P/Halley at minimum phase angle[END_REF] for comet Bowell (1982 I). The reason to adopt the largest value of the linear phase coefficient found in the study of [START_REF] Meech | Observations of Comet P/Halley at minimum phase angle[END_REF] was the large increase of the brightness between Nov 10 and Nov 19, which might be a phase effect (see Table 1). Comparison of the inverse models of the three observations will allow to improve the phase function (see section 4.3.2 and Table 4).

Description of the inversion technique

The direct Monte Carlo calculation provides a set of shells M(x, y; t q , d i ) belonging to the emission time t q and diameter d i and denoting the pΦf contribution of 1 particle, emitted per second. The observed surface brightness b(x, y) expressed in albedo × filling factor units can be described as the sum of the modeled values at that position, multiplied by factors F q,i that are to be de-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
termined by the inversion technique (see eq. 1).

b(x, y)

= Nq q=1 Ns i=1 F q,i M(x, y; t q , d i ) (12)
The coefficients F q,i have the dimension of a particle production rate. In case of a perfect agreement between the direct model and the observation we would have F q,i = Q d (t q )/N s , the total dust production rate of one quasisynchrone divided by the number of shells. Note that in this case all factors F q,1 ...F q,Ns would be equal, as the size distribution within a quasisynchrone has already been taken into account in the particle size distribution function g(t q , d i ). The coefficients F q,i must be ≥ 0, since the dust tail is optically thin. In the inverse model F q,i must be found from fitting the model brightness (right hand side of eq. 12) to the observed brightness (left hand side). From the observed image we may consider only points x, y that contain valid information. We impose two conditions: (i) the selected points exceed a signal-to-noise ratio of three, and (ii) the isophote of this threshold value should not touch the borders of the image. The first condition is essential. The second condition is necessary because close to the edge of the observed field some vignetting may occur.

The typical number of valid pixels N v determined in this way amounts to about thousand (see Table 3 below). In addition the number of participating shells M(x, y; t q , d i ) must be restricted. The factor F q,i of a shell lying totally outside of the valid domain of the observations is undetermined. Such a shell must be removed from eq. ( 12). Shells containing very few valid points could in principle be removed as well but were kept in our algorithm. Shells containing no or only very few valid points most often belong to small particle sizes, as such particles move fast and therefore leave the valid domain comparatively quickly. Shells that contain many valid points but differ very little from each

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
other also cause ambiguity in the inverse solution. Frequently they belong to large particles observed at different times, as such particles move very slowly.

They are not easily recognized and were not removed. After removal of invalid shells the number of valid shells within a quasisynchrone will be less than N s = 51 and we will denote it with N sq ≤ N s . The total number of shells of all quasisynchrones having at least one point within the set of valid pixels

N sh = q N sq ≈ 300 (see Table 3 below).
In the following we define χ 2 relative to to the difference between the number of valid points to be fitted N v and the number of free parameters N f multiplied with σ b (x, y) 2 , the square of the statistical error of the individual image pixels.

For the general inverse problem the number of free parameters N f = N sh , the number of shells. With this definition we can state the inverse problem as the task to minimize χ 2 given by the following equation:

χ 2 = x,y Nq q=1 Nsq i=1 F q,i M(x, y, t q , d i , t obs ) -b(x, y) 2 (N v -N f ) • σ b (x, y) 2 (13)
In this equation the pairs (x, y) are the valid pixels in the observed image. The problem stated by eq. ( 13) for F q,i ≥ 0 is ill posed. Minimization of χ 2 yields unstable results with a large number of coefficients F q,i being negative. Fulle (1989) developed a regularizing formalism which provides stable solutions with only positive values. In the analysis described in this paper we use another approach of solving the ill posed inverse problem. This approach is based on the technique applied by [START_REF] Lemen | Coronal activity in F-, G-, and K-type stars. III -The coronal differential emission measure distribution of Capella, Sigma-squared CrB, and Procyon[END_REF] for derivation of plasma parameters from comparison of observed and theoretical X-ray spectra of the coronae of cool stars. Since we must deal with particle sizes and emission times as two independent quantities we had to generalize Lemen's algorithm to two dimensions.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

The only assumption made in this method is that F q,i can be expressed in the following form:

F q,i = exp K 0 + nq n=1 ns m=1 K n,m C n (Y q ) • C m (X i ) . ( 14
)
where C n (Y q ) and C m (X i ) are Chebyshev polynomials of order n and m, n q is the highest order taken into account for the emission times and n s is the highest order for the particles sizes. The quantities Y q and X i are given by

X i = -1.0 + 2.0 × (i -1)/(N s -1) for 1 ≤ i ≤ N sq Y q = -1.0 + 2.0 × (q -1)/(N q -1) for 1 ≤ q ≤ N q . ( 15
)
i.e. they are distributed in equidistant rows and columns in the square domain If n q and n s are large, eq. ( 14) can approximate any set of positive coefficients F q,i to any desired accuracy. But in most practical cases n q N q and n s N s . The number of degrees of freedom N f equals the number of Chebyshev polynomial coefficients, i. e. N f = n q + n s + 1. As compared to the number of shells N sq it is drastically reduced and eq. ( 14) forces a smooth variation on the coefficients F q,i . If n s = n q = 0 all coefficients will be equal and we arrive at the solution obtained by trial and error with the direct simulation, i.e. the solution having particle size distributions g(t q , d i ) for each quasisynchrone q and a production rate variation ∝ r -2 h . If we insert eq. ( 14) into eq. ( 13) and

-1 ≤ X ≤ 1, -1 ≤ Y ≤ 1 of
A C C E P T E D M A N U S C R I P T
ACCEPTED MANUSCRIPT make use of the expression for the reduced number of degrees of freedom, χ 2 takes the following form:

χ 2 = x,y Nq q=1
Ns,q

i=1 exp K 0 + nq n=1 ns m=1 K n,m C n (Y q ) • C m (X i ) M(x, y, t q , d i ) -b(x, y) 2 (N v -n q -n s -1) • σ b (x, y) 2 (16)
Eq. ( 16) is nonlinear. Therefore the minimization of χ 2 requires an algorithm able to handle non-linear equations. A good initial guess is needed. We constructed the initial guess in two steps. First we weighted the results from the direct Monte Carlo solution with the initial particle size distribution function

g(t q , d i ) ∝ d ν with ν = -3.
3 for all emission times t q used in the simulation of the Nov. 10 and Dec. 03 observation, and ν = -3.2 for all emission times t q used in the simulation of the Nov. 19 image. Secondly we derived a starting value for K 0 . The latter was performed by requiring the total brightness in the model to be equal to the total brightness in the observation, and by setting all other K n,m (n, m > 0) = 0.

With this initial guess the solution with order (1, 1), i. e. order 1 in time and order 1 in size was calculated. Higher orders were derived recursively, i.e.

the solution having one order less either in time or in size was used as initial guess for the higher order. At orders higher than three the solutions become unstable.

There is no doubt that eq. ( 13) does not have a unique solution, i.e. the inverse problem cannot be solved uniquely. In contrast to the method used in the past by Fulle and his coworkers, the new method provides several solutions depending on the applied order. This allows to investigate how sensitive the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
solution selected by the method is to the choice of the order. We cannot get an overview over the full set of available solutions, but at least we can investigate a small subset of them. When we explicitly write down the low order terms of eq. 14 we notice that, because of the presence of the exponential function, the logarithm of the function to be fitted will be a polynomial of the coefficients to be determined by the method. Let us first consider particle size. In the direct method the size shells are logarithmically spaced but the corresponding coefficients enter eq. ( 15) linearly. Therefore, if we plot the logarithm of particle number per size bin versus the logarithm of size, in first order we will get straight lines, i.e. power law size distributions. In second order the logarithm of particle number will be a parabolic function of the logarithm of size. Similarly a first (second) order in time will create solutions with linear (parabolic) deviation of the logarithm of the production rate from the inverse square law in heliocentric distance applied in the direct simulation. This will become more clear in the next section when we will discuss the application of our inverse technique to the three images of comet WM1.

Application of the model

The calculated inverse solutions and their χ 2 values

In order to interpret the three observations of Nov 10, Nov 19, and Dec 03 we applied the above described method of minimization to eq. ( 16) using the shells M(x, y, t q , d i ) calculated from the three images as input. The minimization of χ 2 was performed with the aid of the IDL realization of Marquardt's gradientexpansion algorithm. We determined the orders 1 to 3 in n q (time) and 1 to 2 in n s (size) of the polynomial expansion ( 14). The achieved χ 2 values for the was made to find a single solution that fits the three images simultaneously.

Determination of empirical values of the brightness phase curve

Minimization of eq. ( 16) with the method described in section 4.2 allows us to determine the coefficients F q,i that, according to eq. ( 12), link the shells M(x, y; t q , d i ) with the observed surface brightness b(x, y). But the shells M do not only depend on the parameters of the direct model but also on the albedo phase function product pΦ(φ) (eq. 11), which has not been specified yet. As described in section 4.1.3 we adopted for the dust grains of comet WM 1 a geometric albedo of 4%, the albedo of the nucleus of comet Halley. The phase function was written in the form [START_REF] Meech | Observations of Comet P/Halley at minimum phase angle[END_REF]) Φ(φ) = 10 -0.04 γ φ , with the linear phase coefficient γ = 0.035. Table 4 provides in its first row the phase angles at the three nights of observation and in its second row the product albedo-phase function as derived from this formula. Now, with solutions of the inverse problem at hand, we can calculate for each of the three observations "particle area production rates" using the adopted values of the albedo-phase function product (Table 4, second row) and check, if the area production rates derived from different observations agree for the quasisynchrones common to more than one observation. We found that they were systematically different. Therefore we readjusted the albedo-phase function product of the three observations in order to improve the agreement as much as possible. The readjusted values are given in the third row of Table 4. The (2) a surge of the phase function at about 12 • degree phase angle, the minimum value observed. Theoretical phase function models do not predict such a high brightness surge at a phase angle as large as 12 • . We note that the adjustment of the phase function automatically takes care of the fact that our three observations were taken through three different continuum filters (Table 1). The wavelength dependence of albedo enters the values of the albedophase function product presented in Table 4 but cannot be separated from it. Therefore we should look at the adjusted values as parameters needed to compare the results derived from the different observations, but should not interpret them in terms of scattering properties of the cometary dust grains. 3. Here and in the following the notation (i, k) indicates order i in particle emission time and order k in particle size. All combinations provide a satisfactory agreement. If the emission time order increases from 1

Results of the inverse model and discussion

Comparison of observed and modelled images

to 2 the agreement of observation and model is significantly improved. This is, however, not the case if the emission time order rises from 2 to 3. The

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

in situ measurements of the dust particle distribution in comet 1P/Halley have revealed power law size distributions in wide size ranges [START_REF] Mazets | Comet Halley dust environment from SP-2 detector measurements[END_REF][START_REF] Mcdonnell | The dust distribution within the inner coma of comet P/Halley 1982i -Encounter by Giotto's impact detectors[END_REF]. Therefore particle size distributions significantly different from a power law like those generated by particle size orders > 1 do not appear realistic (Fulle, 2004). Nevertheless, besides of the solutions with particle size order = 1 we will also present the inverse solution of orders (2,2)

as an example of a size distribution generated by a second order polynomial of log(d).

Because of the great importance of particle size distributions in the form of a power law, in the next section we recall some properties of the integrals over power law distributions that are important for the interpretation of our results.

Properties of integrals over power law distributions

In first order of particle size our method yields power law particle size distributions (section 4.2). Besides of the production rate of the number of particles per second, i. e. the direct integral over the particle size distribution, two other production rates are of importance, the dust mass production rate and the "particle area production rate" we used in section 4.3.2 to derive an empirical phase law (see Figure 12). In the particle size range investigated in this paper the scattering cross-section of the particles is very nearly proportional to their geometrical cross-section (section 3). Therefore we can write these three types of production rates as integral over a power law distribution:

Q n = c n (d 0 ) • dmax d min d n+α d d, n = 0, 2, 3, (17)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

where Q 0 , Q 2 , and Q 3 are the production rates of number, area and mass, respectively. c n (d 0 ) is a reference value of differential particle number, area or mass per cm and s at the reference diameter d 0 . d max and d min denote the maximum and minimum cutoff diameter of the size distribution. α is the power index. Note that in this paper, like in Fulle (1999Fulle ( , 2004)), the definition of the power index includes its sign, i.e. α usually is negative. The integral can be easily solved in terms of analytic functions. We obtain:

Q n = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ c n (d max ) • d max • 1 n+1+α 1 -d min dmax n+1+α : α > -n -1, c n (d min ) • d min • -1 n+1+α 1 -d min dmax -n-1-α : α < -n -1, ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ n = 0, 2, 3.( 18
)
In this equation the reference value c n defined above is taken at the maximum and minimum cutoff diameters d max and d min . We have introduced a case distinction into eq. ( 18) to assure that the power associated with the ratio Then the minimum diameter is much higher: d min = 100.0 μm. In any case we have d min /d max ≤ 0.2. If we consider 0.2 1.5 = 0.09 negligible as compared to unity we can derive the following conclusion from eq. 18:

For

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ α > -n -1 + 1.5 Q n depends only on the reference value at d max , α < -n -1 -1.5 Q n depends only on the reference value at d min , ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ n = 0, 2, 3.(19)
In the bottom row of Figure 14, for the order combinations that have particle size order = 1, i.e. power law size distributions, the values of the power index of the particle size distribution derived from the inverse model are displayed The number and mass production rates will be discussed in section 5.4.

Derived particle size distributions

The particle size distributions and their temporal evolution is shown in Figure 13. The three panels refer to the three observing days. For each observing day the size distributions of the individual quasisynchrones are stacked on top of each other. The individual quasisynchrones cover different size ranges. Small particles have higher initial velocities and are more strongly accelerated by solar radiation pressure. Therefore early quasisynchrones lack small particles, if the time between their emission and the observation was long enough to allow them to leave the field of view (see also Figures 678). E. g., particles with sizes smaller than about 100 μm emitted more than 228.4 days before perihelion (quasisynchrone 1, Table 2) have left the field of view of the Nov 10 observation.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

The particle size distributions derived from the Nov 10 image do not vary strongly, neither with time (quasisynchrone number) nor with the particle emission time order. For order combination (1,1) the particle size distributions are very nearly parallel. With respect to the reference points on the right side the curves move downward with time (increasing quasisynchrone number), i.e.

for all particle sizes the number of particles decreases with time. This is at variance with the fact that the comet approaches the Sun and therefore particle production is expected to increase. We will encounter the same effect in the solution derived for Dec 03. As will be explained in the discussion of Dec 03 it is very likely caused by overlap of shells of heavy particles released shortly before the observation. The other order combinations (2,1), (3,1) and (2,2) agree more closely among each other than with (1,1). The early quasisynchrones show a slightly steeper particle size distribution and an increase of the particle number in the distributions. The later quasisynchrones (no. 5 -8) are flatter and the particle content decreases as it does for order combination (1,1).

On Nov 19, because of the anti-tail geometry with the Earth close to the cometary orbit plane, the shells are not located as favorably as on Nov 10.

At the large particle end of the particle size distribution, for all order combinations including (2,2), the number of particles rises with time. But there are large instabilities in the power index of the particle size distribution for order combinations (1,1) to (3,1) in the early and late quasisynchrones. Consider quasisynchrones 8 and 9 for order combinations (2,1) and (3,1). The power indices are large and, according to our theorem (eq. 19, section 5.2) the brightness contribution comes mainly from the large particles. Therefore the size distributions of both orders agree for the large particles but differ (are unstable) for the small particles. Conversely, for quasisynchrone 2 and order

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

combinations (1,1) and (3,1), the power indices are small and, according to our theorem, small particles contribute most to the brightness. Consequently the particle size distributions of the two orders meet at the small particle end of the size distribution and differ (are unstable) for the large particles. To a lesser extent this applies to quasisynchrone 1 as well. The power index of order combination (3,1) is -12.2 (out of the plotting range). Evidently the instability is so large here that the two order combinations cannot anymore be compared directly and a small difference in particle numbers remains at the small size end of the size distribution.

As compared to the November observations, on Dec 03 our view on the comet has significantly changed. The phase angle has increased from 12.5 • to 44.6 • (Table 1). While in November the viewing direction was nearly parallel to the tail axis, on Dec 03 we see the comet much more "from the side". Moreover, the comet has approached the Earth. It reached its minimum distance of 0.3164 AU one day before the observation. Because of this, as compared to the November observations, only a rather small part of the comet is contained in the field of view. Most of the shells of the early quasisynchrones are outside of the field of view. Therefore only quasisynchrones 6-11 could be used for the inversion procedure.

On Dec 03 for order combinations (1,1), 2,1) and (3,1) there is a trend for the power index to decrease with time (Figure 14). Consequently, in accordance with our theorem, if we compare different order combinations and quasisynchrone numbers in Figure 13 we notice a gradual shift. For quasisynchrone 6 the number of large particles does not change with order but, because small particles do not contribute to the brightness, the number of small particles in the distribution is unstable. For quasisynchrones 11 and, to a lesser extent,

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
quasisynchrone 10, the small particles contribute most to the brightness and therefore the particle size distribution is stable for the small particles and unstable for the large ones. Order combination (2,2) also follows the described trend.

Probably the increase of small particles in the size distribution is at least partly real. The shape of the observed dust tail as compared with the pattern of synchrones and syndynes (Fig. 3) suggests already the high abundance of small particles: The dust tail bends down but the bend is not strong enough to favor the central syndyne with β = 0.01. Instead syndynes of higher β and smaller particle size dominate the picture. But the shells of recently emitted large particles may as well not be selected for reasons inherent in the inversion method (see Figure 8). First of all the shells of recently emitted large particles overlap strongly. Second, they are sharply peaked. It is likely that the peaks are more pronounced than the near-nuclear intensity increase in the observed image, as the observed peak is smeared by seeing effects and guiding errors.

Last not least, the shells in question are located at the border of the intervals of both emission times and particle sizes where our inversion method is least accurate. We will return to this point in the next section.

Production rate of particle number and mass: Comparison of different observation days and order combinations

The temporal evolution of the production rates of particle number and particle mass are important quantities. Having at hand the particle size distributions we can calculate the temporal evolution of the production rates of total particle number and mass for the different order combinations and observing days.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

The resulting production rates are shown in Figure 14. The bottom panels of this figure display the power index of the particle size distribution discussed already in sections 5.2 and 5.3. The upper and center panels show particle and mass production rates (a density ρ = 1.0 g cm -3 has been adopted).

As was said above, the smaller particles are emitted from the nucleus with higher initial velocity and they are more strongly accelerated by radiation pressure. Therefore small particles are missing in the older quasisynchrones.

To compensate for this time dependence of the particle size range contained in the observed images we have calculated two production rates. Minimum production rates were derived for the particle size range of large-sized particles present in all quasisynchrones (100.0 ≤ d ≤ 524.8μm, d is particle diameter).

In addition a maximum rate was calculated by extrapolating the particle size distributions to the maximum range of particle sizes (4.786 ≤ d ≤ 524.8 μm)

using the available functional dependence of the size distribution. The different observing days are plotted with different line styles and symbols. The plotting style is the same for the minimum and maximum rates.

For the different observing days the particle number and mass production rates should be the same but they agree at best by order of magnitude. We must remember here that, because of the poorly known phase curve and spectrophotometric gradient (see Section 4.3.2), in Figure 14 (upper and center panels) we can shift the curves of the different observing days up and down with respect to the curves of the other two observations by some amount.

But the differences in Figure 14 concern the shape of the curves rather than their relative location. Moreover, contrary to expectation (see Section 5.3), the dust production rates do not always increase with decreasing heliocentric distance. While the order combinations (2,1), (3,1) and (2,2) produce a

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
similar behaviour of the production rates the order combination (1,1) differs significantly.

To understand the behavior of the maximum and minimum production rates in Figure 14, eq. ( 19) derived in section 5.2 is again useful. We start with the production rate of particle number and consider again the case of Nov 19. As the lower limit of the particle size distribution, because of particles leaving the field of view of the observation, varies between 4.786 and 100.0 μm, a reasonable particle number production rate can only be obtained if this production rate is determined by the large particle cutoff size of the distribution.

According to eq. ( 19) this is the case if the power index α > +0.5. For order combination (2,1) the maximum power index is reached at the earliest and latest quasisynchrone, where it amounts to -0.18 and -1.05, respectively, still quite close to the critical value α = -1. Nevertheless the maximum and minimum particle production rates for quasisynchrones 1 and 9 are close to each other. For the other quasisynchrones in order combination (2,1), and for many quasisynchrones of other observation days and order combinations as well, the power indices α ≤ -3 and the minimum and maximum particle number production rates differ by many orders of magnitude. Therefore, because of the variable lower cutoff limit of particle size, the errors in particle number production rates can reach several orders of magnitude, and we just cannot determine them with our models. Fortunately, the situation is not quite as bad when we consider the much more important dust mass production rates.

For the dust mass production rates the critical power index α = -4. Most of the power indices in Figure 14 are larger than that but not many exceed the limit of -2.5 of eq. ( 19). As can be seen from the center panels of Figure 14 the minimum and maximum mass production rates in order combinations (1,1)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
and (2,1) rarely differ by more than a factor of three. But there are still much larger differences between the different observations. Particularly worrisome is the strong decrease of mass production toward smaller heliocentric distances seen in Figure 14 for the observing days of Nov 10 and Dec 03. Is this effect real?

A 'Hearn et al. (1995) have derived dust production rates for a large number of comets using aperture photometry of Af ρ. They find for dynamically new comets before perihelion d log(Af ρ)/d log(r hel ) = +0.03, i.e. a very small decrease of dust production rate with decreasing heliocentric distance, somehow in agreement with our result. We must note, however, that A'Hearn et al.

(1995) simply take the measured Af ρ values and do not consider the dependence of emission velocity on heliocentric distance. In our paper we have used empirical values of the particle ejection velocity derived from the observations (Table 2, 10th column) or the power exponent of -1.5 (eq. 9). If we apply this power law to the exponent +0.03 [START_REF] A'hearn | The ensemble properties of comets: Results from narrowband photometry of 85 comets, 1976-1992[END_REF] we get a power of -1.47, i.e. a rather strong increase of the dust production rate with decreasing heliocentric distance. We must therefore conclude that from Af ρ aperture photometry we must expect a significant increase of dust production with decreasing heliocentric distance even for new comets. This conclusion is in agreement with the discussion about Af ρ by Fulle (2000Fulle ( , 2004)).

If we compare the results obtained with emission time order 1 with those of higher emission time order we notice that in the Nov 10 observation higher emission time orders provide an increase of the dust production rate with decreasing heliocentric distance at least at larger heliocentric distances. Evidently, as already mentioned in Section 5.3, an emission time order of 1 does not yield an appropriate inverse solution for Nov 10. Apart from the fact that

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
the agreement between model and observation (Figure 9) is somewhat better for emission time orders > 1 this is also indicated by the closer similarity of the total mass production rates of Nov 10 and Nov 19 derived for emission time orders > 1 in Figure 14. Fulle (2000) has constructed inverse numerical models of comet 46P/Wirtanen for the time before its perihelion in March 1997. Most of his models indicate that during approach to perihelion the dust mass loss rate decreased by about a factor of 2. After careful consideration of possible errors he has come to the conclusion that this decrease is real and accompanied by a decrease of the dust size distribution power index to a value < -4. The case of Dec 03 is similar to his one but, even if we consider only emission time order = 1, much more dramatic. With a power index < -4 the mass is determined by the lower cutoff size of the distribution. Therefore, to rise the mass loss we must decrease the lower limit of the size distribution and include smaller particle sizes in our inverse model. From Figure 3 we cannot exclude the presence of particles of diameters smaller than 4.786 μm, in particular as the syndynes of small particles accumulate in the anti-solar direction. But our attempts to produce inverse solutions with a larger amount of small particles were not successful.

The agreement between observation and model was always worse than for the solutions provided here. All model solutions with a smaller lower limit of the particle size resulted in the appearance of a dust tail pointing exactly in antisolar direction. This is not observed in our images. Therefore it is doubtful if inclusion of smaller particles can produce a dust mass loss that increases with decreasing heliocentric distance and at the same time still provides a good agreement between observation and model. Another possibility to increase the production rates at times close to observation time was mentioned already in

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
Section 5.3. There we gave reasons why our inversion method may not select the most recently emitted shells of large particle size. Replacing a shell of small particles with an equally bright shell of larger particles will increase mass content without increasing the scattering area, i.e. this is exactly what is needed to improve the mass balance.

Not only on Dec 03 but also on Nov 10 we have the phenomenon of a decrease of the mass production rate with decreasing heliocentric distance. But at the time when the production rate derived from the Nov 10 image decreases, the rate derived from the Nov 19 image increases. This instability adds to our suspicion that the decrease may not be real.

In view of these arguments we must conclude that the mass production rates and even more the production rates of particle number shown in Figure 14 must be interpreted with great caution, in particular if small particles are so frequent that they carry most of the particle mass.

Assessment of the new inversion method

After the pioneering work by Finson and Probstein (1968b) Fulle (1987, 1989) was the first to construct a dynamical dust tail model adequate to determine properties of real comets. According to Fulle (1999Fulle ( , 2004)), but expressed in terms of our paper, his model (1) computes Keplerian orbits by sampling millions of dust grains, (2) uses shells described by sizes and emission times,

(3) takes into account anisotropic dust ejection, (4) provides for each dust size a dust ejection velocity from the coma that is a mean value of a wide velocity distribution, and (5) avoids the trial-and-error procedure typical of the original

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
Finson-Probstein model by means of an inverse ill-posed problem theory. In this section we compare Fulle's approach with ours, addressing points (4) and

(5).

The treatment of the ejection velocity of dust particles in our model is very similar to the one used by Fulle. The velocity is dependent on particle size, time of ejection, and emission direction (anisotropy). Like Fulle we assume separability of variables, i.e. the dust grain velocity is a product of three functions depending on time (or heliocentric distance), particle size and emission direction, respectively (eq. ( 6)). In our model (Section 4.1.1) the mean dependence of the velocity on heliocentric distance is well described by a power law with index = -1.5. This value is between the values derived from modeling the coma and tail of comet 46P/Wirtanen (Fulle, 2000) and from parametrization of the results of hydrodynamical coma models by Crifo and Rodionov (1997).

But, like Fulle, in addition we derive the time dependence of the velocity from the data images using trial and error. We describe the dependence of particle emission velocity on particle size by a power law for particles larger than a given size, with a smooth transition to a constant limiting value for the smaller grains (eq. ( 8)). But the observations considered in this study are well described by relatively large particles, so we are working in the power law range with an index = -0.5. Other models (Fulle, 1999(Fulle, , 2004;;[START_REF] Kimura | On the structure of cometary dust tails. Chinese A C C E P T E D M A N U S C R I[END_REF][START_REF] Waniak | A Monte Carlo approach to the analysis of the dust tail of Comet P/Halley[END_REF][START_REF] Waniak | Nuclear dust emission pattern of Comet Wilson 1987VII[END_REF] derive larger indices like -0.25. In this case the velocity of large grains relative to that of small ones is higher than predicted by one-dimensional coma models of the dust-gas interaction and perhaps better described by our eq. ( 8). Concerning the emission directions Fulle et al. (1992) assume a cone with halfwidth w (0 < w < π), pointing toward the Sun. In the cone the velocity function is unity and zero elsewhere. In our model the et al., 1998) or if the comet behaves erratically (Bonev et al., 2002). To solve the inverse problem Fulle (1987Fulle ( , 1989) ) used the method described by [START_REF] Tikhonov | Sublimating components in the coma of comet C/2000 WM 1 (LINEAR)[END_REF] and imposed as regularization condition that the particle size distribution be slowly time variable and the dust production rate have an inverse square dependence on heliocentric distance. Fulle (1987) did not want to use the condition that the coefficients F q,i (eq. 12) should be positive.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
According to Fulle (1987), if too high ejection velocities are adopted, the dust shells will overlap each other excessively and the coefficients F q,i will sometimes be negative to balance the excess of overlapping. In our view overlap of shells is unavoidable. Even though a dust tail is to some extent confined to the orbit plane of the comet the problem is three-dimensional and the shells overlap in areas where the motion of the particles is close to the line of sight. Old, slowly moving shells of large grains can overlap with more recently released shells of smaller grains.

Besides of the Tikhonov method other methods to solve the inverse problem have been applied. [START_REF] Moreno | Dust in Comet 67P/Churyumov-Gerasimenko[END_REF] employ the singular value decomposition technique (see e.g. Press et al., 1992, , p. 402). In their dynamical model determining the three-dimensional outflow from the nucleus [START_REF] Waniak | Dust Emission for Comets Shoemaker-Levy 1991a1 and McNaught-Russell 1993v[END_REF] use the Richardson-Lucy (Richardson, 1972;[START_REF] Lucy | An iterative technique for the rectification of observed distributions[END_REF] algorithm that helped to sharpen the early Hubble Space Telescope images. This algorithm is designed to ensure that the coefficients F q,i will be positive. Its merit lies in the fact that it is derived from maximum entropy considerations that are appropriate to our problem as well.

Our mechanism [START_REF] Lemen | Coronal activity in F-, G-, and K-type stars. III -The coronal differential emission measure distribution of Capella, Sigma-squared CrB, and Procyon[END_REF]) is conceptually simple. Despite its nonlinearity, solutions could be found by a straightforward application of Marquard's gradient-expansion algorithm (see Section 4.2). We can select different orders

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
of the polynomials describing the dependence of the solutions on particle size and temporal variation of production rate. As an advantage, the method allows to restrict the inverse solution to power law size distributions from the very beginning by selecting order 1 in particle distribution. But there are disadvantages. They are related to the inherent instability of the representation of the flux coefficients F q,i by Chebyshev polynomials (eq. 14), in particular in the two-dimensional case. When higher orders are selected, there is a tendency that eq. ( 14) favors a few coefficients F q,i at the corners of the coefficient domain of eq. ( 15) and connects them to a few Chebyshev coefficients C n and C m . Our method, in particular when used at higher orders, favors "island" solutions where the observations are reproduced by massive contributions from a few selected shells. One example we have encountered in section 5.3: Quasisynchrone 1 of Nov 19 is attributed a very low power index of -12.2 in order combination (3,1) and a very high power index of (-0.18) in order combination

(2,1). I. e., when going from order (3,1) to order (2,1) a single shell with very small particles is replaced by a single shell of very large particles, of course with rearrangment of the shells of intermediate sizes in this and in other quasisynchrones. Contrary to that, regularization methods based on maximum entropy probably favor solutions were many shells participate. But the fit to the observation is similarly good for all these solutions! We should not forget, however, that in this paper we deal with preperihelion observations where the overlapping of shells is much more severe than after perihelion. Also observations in a large field of view are better suited for dynamical models than narrow-field observations. In our case this applies particularly to the Dec 03 observation.

From a practical point of view we give the advice not to use Lemen's method

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
for orders ≥ 2 in particle size and ≥ 3 in emission time. As suggested by Waniak (2007, private communication), one could perhaps apply the mechanism iteratively in low orders. The result of a first iteration could be implemented into the functions g(t q , d i ) of eq. ( 11) of the direct solution and

Lemen's method could be used a second time.

Changes in the size distribution of the dust particles released by comet WM 1 , comparison with other comets

As apparent from Section 5.3, as comet WM 1 approached the Sun, between Nov 19 and Dec 03 (1.423 > r h > 1.186) an increasing number of small particles appeared in the dust coma and tail of WM 1 . In this Section we compare comet WM 1 with comets C/1999 S4 (LINEAR) and C/1996 B2 (Hyakutake) which were observed before perihelion and in which transitions in the particle size distribution to an increased number of small particles were also observed.

Comets WM 1 , C/1999 S4, and C/1996 B2 certainly belong to the class of "nearly isotropic comets" as defined by [START_REF] Levison | Comet Taxonomy[END_REF]. The question if they were observed entering the inner solar system for the first time cannot easily be answered, however. For a number of nearly isotropic comets with well-known orbits Dybczyński ( 2006 Comet C/1999 S4 dissolved in August 2000 [START_REF] Boehnhardt | Comet splitting -Observations and model scenarios[END_REF][START_REF] Boehnhardt | Split Comets[END_REF]) at a heliocentric distance very close to its perihelion distance of 0.77 AU. During the days before disintegration Bonev et al. (2002) observed two strong brightness outbursts that were associated with the emission of small dust grains. Bonev et al. (2002) concluded that "in particular in the second, final outburst the excess small particles could originate from fragmentation of fresh larger particles". After disintegration, however, a "sandbank" of large particles remained indicating that not all dust was broken up into small grains. The two outbursts were associated with peaks of the water production rate (Figure 2 of [START_REF] Mäkinen | Water Production of Comet C/1999 S4 (LINEAR) Observed with the SWAN Instrument[END_REF]. During the final break-up Bockelée-Morvan et al. a peak in the water production rate observed on April 14.6 [START_REF] Hicks | Spectrophotometry and the Development of Emissions for C/1996 B2 (Comet Hyakutake)[END_REF][START_REF] Biver | Spectroscopic Monitoring of Comet C/1996 B2 (Hyakutake) with the JCMT and IRAM Radio Telescopes[END_REF][START_REF] Combi | Temporal deconvolution of the hydrogen coma[END_REF]. It must be distinguished from other outbursts observed in this comet on April 9 and around perigee (March 25, 1996), for which no information about the dust particle size distribution is available.

Was the change in the particle size distribution found in comet WM 1 associated with other forms of cometary activity? The possible existence of two sublimating components in the coma of this comet, derived from observations of the albedo filling factor product Af ρ on Dec 02 and Dec 04 2001 by Tozzi et al. ( 2004), can perhaps be alternatively explained by the change of particle size distribution resulting from our dynamical modelling. Figure 15 shows brightness data (courtesy of International Comet Quarterly), and HCN [START_REF] Biver | Radio wavelength molecular observations of comets C/1999 T1 (McNaught-Hartley), C/2001 A2 (LINEAR), C/2000 WM 1 (LINEAR) and 153[END_REF] and CN [START_REF] Schleicher | Photometry and Imaging of Comet C/2000 Wm1 (Linear)[END_REF][START_REF] Lara | Gas and dust in Comet C/2000 WM1 during its closest approach to Earth: Optical imaging and long-slit spectroscopy[END_REF] production rates of comet WM As small grains have a larger area as compared to large grains of the same total mass and scatter the solar light more efficiently, the brightness increase probably is related to the change in particle size distribution. The CN data [START_REF] Schleicher | Photometry and Imaging of Comet C/2000 Wm1 (Linear)[END_REF] are smooth but have low temporal resolution. The HCN data [START_REF] Biver | Radio wavelength molecular observations of comets C/1999 T1 (McNaught-Hartley), C/2001 A2 (LINEAR), C/2000 WM 1 (LINEAR) and 153[END_REF] cover the period from Nov 23 to Dec 08 very close to our period of interest and also indicate a smooth behavior of the outgassing rate. Therefore, as opposed to the events observed in comets C/1996 B2 (Hyakutake) and C/1999 S4 (LINEAR), in WM 1 no increase of the gas

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
production rate accompanied the appearance of small particles. Another, very large increase of brightness occurred shortly after perihelion, when the comet was at a heliocentrc distance of 0.55 AU, i. e. at the distance of the Hyakutake event. But it is not known if it was associated with an increased release of small dust grains. Published observations of comet WM 1 in the thermal infrared spectral region are not available in the literature.

In-situ observations of cometary dust indicate that dust fragmentation is continuously going on in comets 1P/Halley [START_REF] Keller | Hydrodynamic implications of particle fragmentation near cometary nuclei[END_REF][START_REF] Oberc | Disintegration of Dust Aggregates as Origin of the Boundaries in Halley's Coma: Derivation of the Sublimation Parameters[END_REF] and 81P/Wild 2 [START_REF] Clark | Release and fragmentation of aggregates to produce heterogeneous, lumpy coma streams[END_REF]. [START_REF] Oberc | Small-scale dust structures in Halley's coma II. Disintegration of large dust bodies[END_REF] has given a comprehensive treatment of dust disintegration as observed in comet Halley: Particles with a mass close to the maximum mass that can be lifted from the nucleus fragment in a multi-step process. The particles are of fractal nature. The fragmentation mechanism is sintering combined with sublimation of an organic glue with a latent heat of ≈ 80 kJ mol -1 , corresponding perhaps to cometary polymers [START_REF] Boehnhardt | The possible role of organic polymers in the structure and fragmentation of dust in the coma of Comet P/Halley[END_REF][START_REF] Huebner | Polyoxymethylene in Comet Halley[END_REF]. Sublimation of ices (water ice or more volatile ices) will lead to catastrophic disruption of the grains and is therefore too fast for the fragmentation observed in comet Halley. It may, however, be the mechanism of choice to explain the events in comets C/1996 B2 (Hyakutake) and C/1999 S4 (LINEAR) discussed above when partial break-up of a nucleus releases volatile matter from its interior.

We conclude that the appearance of small grains in comet WM1 during its preperihelion approach to the Sun (1.423 > r h > 1.186) was not accompanied by a major gas release. Attempts to implement small fractionated grains with low emission speeds in addition to those described by equation ( 8) did not significally improve the agreement between model and observations. Therefore we cannot be sure from our work if the small grains were released directly vector. The method works with Chebyshev polynomials of different orders in emission time and particle size. The lower the order the stricter is the enforcing of the smoothness. In lowest order of particle size the method yields power law size distributions. Orders are written in the form (n, m), where n is the order of emission time and m the order of particle size. Models with 1 ≤ n ≤ 3 and 1 ≤ m ≤ 2 were calculated, of which order combinations (1, 1), (2, 1),

(3, 1) and (2, 2) were investigated in more detail. The following results were obtained.

(1) The three observations are well fitted by the models (Table 3,Figures [9][10][11]. The goodness of fit increases somewhat as the order n increases from 1 to 2. For n ≥ 2 all combinations fit equally well.

(2) But the particle size distributions and their time dependence disagree significantly in models belonging to different orders (Figure 13). For solutions with m = 1 (power law size distributions) this implies that the power indices sometimes vary strongly for different values of the emission time order n (Figure 14 bottom row). All these differences are the natural consequence of the non-uniqueness of the solution of the inverse problem. They are particularly large at the earliest and latest emission times considered for an observation. This can be attributed to the enforced smoothness requirement in time that is more restrictive at the center of the interval than on its start and end points. In the center of the emission time ranges of the three observations the power indices α are in the range of accepted power indices -4 < α < -2 (Jockers, 1997b;Fulle, 1999Fulle, , 2004) ) for all investigated emission time orders n. At the end points, however, the power index α takes values -12.2 < α < 0.51. Fulle (2000)). This decrease can hardly be real. In our opinion it is caused by the fact mentioned already in the previous item: our method may not properly take into account particle shells of large particles released shortly before the observation.

(5) Dust mass production rates are more accurately predicted by dynamical models than particle number production rates (Section 5.4), because number production rates are more strongly affected by loss of small fast moving particles than mass production rates.

(6) Our inversion method [START_REF] Lemen | Coronal activity in F-, G-, and K-type stars. III -The coronal differential emission measure distribution of Capella, Sigma-squared CrB, and Procyon[END_REF]) is conceptually simple and, despite of its nonlinearity, not difficult to apply. If the particle size order m = 1 is selected, power law size distributions are derived directly. velocities derived by [START_REF] Crifo | Hydrodynamic models of the collisional coma[END_REF] and by [START_REF] Gombosi | A heuristic model of the Comet Halley dust size distribution[END_REF], respectively. Thin lines plotted with different line styles describe velocities calculated for different gas production rates and dust/gas production ratios (eq. 8). The thick line shows the ratio of radiative to gravitational force β (eqs. 4 and 5).

  time of observation. As b and M are images with the number of image points (pixels) exceeding the dimension of the flux vector, eq. (1) can in principle be inverted. But the inversion is not straightforward. Because the number of dust particles is positive the elements of the flux vector F must be positive or zero. Because of the finite ejection speed of the dust grains from the nucleus and sometimes also because of the unfavorable location of the cometary orbit with respect to the observer, particles of different sizes are present along a single line of sight through the cometary dust cloud, i.e. the individual shells forming the model image vector M overlap. This makes the solution of eq. (

  Instead, authors have solved eq. (

  2000 WM 1 (LINEAR) was observed with the 2m RCC (Ritchey-Crétien-Coudé) telescope on Pik Terskol from Nov 10 to Dec 6. In this paper we apply our dynamical model separately to three continuum images of comet WM 1 , obtained on Nov 10, Nov 19 (with anti-tail) and Dec 3, 2001 as part of this campaign. All three images contain information about the dust grains released before Nov 10, i.e. before the first observation. Therefore the availability of models of three separate observations spaced by 9 and 14 days provides an opportunity to investigate the accuracy of the models and, in particular, of the inversion method.

  Section 4.1 is devoted to what we call the direct model, i.e. to a description how we determine the initial velocities and how we calculate the shells, i. e. the model image vector M. The inversion technique is described in Section 4.2 and its application in Section 4.3. Section 5 deals with the results (Figures 9 -14) and their discussion. A particular issue addressed in several sections is the question of particle albedo and of the phase function of brightness (Sections 4.1.3 and 4.3.2, Figure 12, Table

  Fig. 2. Note that the part of the orbit already passed by the comet is projected toward the Sun. Most of the projected dust tail (see Fig. 1 below) belongs to the anti-tail. Only the young tail is projected on the anti-solar side, but in three dimensions it is pointing nearly away from the observer, i. e. the line of sight passes through an extended part of the dust tail and particles of many sizes appear superimposed in the comet tail image.

Figure 3

 3 Figure 3 presents selected syndynes and synchrones superimposed on contour maps of the images obtained on Nov 10, Nov 19, and Dec 03, respectively. The Sun is to the left (negative x-axis). A thin straight line marks the cometary orbit and an arrow the direction of cometary motion. The contours are labeled with the decadic logarithm of the pΦf values. The β ratios of the six drawn syndynes (dotted) are provided in the figure heading. They differ by a factor 10 0.4 = 2.51. The syndyne with the highest β value (smallest particle diameter)

  time steps comprising the quasisynchrone are equidistant in true anomaly and also the time intervals of the quasisynchrones themselves are (with small exceptions) equidistant in true anomaly instead of time. Like real synchrones also the quasisynchrones consist of particles with different sizes (β values). In our model a quasisynchrone contains n s = 51 shells belonging to 51 β values, logarithmically distributed in the range 0.0016 ≤ β ≤ 0.016 (see

  10 4 model particles, independent of the particle diameter d i of that shell. After ejection the particles move along their orbits until the time t obs , when the particles are "observed". To determine the "observed" surface brightness of the shell we must calculate the position of the model particles at time t obs and collect them into model "detector bins" projected into the tangential plane passing through the cometary nucleus [x, x+Δ x] and [y, y+Δ y].

  figure headings. These figures are discussed in more detail in section 5.3.

  the products of Chebyshev polynomials in X and Y . As some shells with no overlap to the valid image domain are removed, the domain -1 ≤ X ≤ 1, -1 ≤ Y ≤ 1 is not evenly filled with anchor points for the coefficients. It is more important that the coefficients of the same shell in different quasisynchrones are on the same straight line in the domain of Chebyshev polynomials.

  the different orders are listed in Table3together with other relevant information. Inclusion of larger orders of Chebyshev polynomials lead only to a very small further decrease of χ 2 . At the same time it introduced oscillations in the final solution that have no physical meaning. No attempt

  rates are shown in Figure 12 for the case of orders one in particle emission time and size distribution. The adjusted values imply (1) no change of the phase function between 45 • and 18 • phase angle and

Figures 9 ,

 9 Figures 9, 10, and 11 show the comparison of the observations (full lines) and the model (dashed lines) for the three nights of Nov 10, Nov 19, and Dec 03, respectively. The six panels of each figure correspond to the order combinations of Table3. Here and in the following the notation (i, k) indicates

d

  min /d max is always positive. The maximum cutoff diameter used in this paper d max = 524.8 μm and the minimum diameter d min = 4.786 μm. In some old quasisynchrones the small particles have left the field of view of the image.

  quasisynchrones as function of heliocentric distance. Let us consider the area production rate. For a power index -4.5 < α < 1.5 individual shells of all particle sizes contribute significantly to the observed brightness. If in Figure14power indices from different observing days at similar heliocentric distances are close to each other, they are in this range. We conclude that in this range of power indices our inversion method is reliable. For a power index outside this range the brightness is determined either by the large or by the small particles. It is obvious that no information can be derived on particles that do not contribute to the brightness of an observed image. The results on such particles are then only based on the smoothness requirement imposed to regularize the inverse problem.

  ) has calculated not only original aphelia but also perihelia of their previous returns, taking into account the gravitational influence of the galactic disk and of nearby perturbing stars. None of our three comets appears in his work. Comet C/1999 S4 has an original inverse semimajor axis 1/a original = +0.000024 (MPEC 2000-N15) and therefore can be considered dynamically new. From high-dispersion spectroscopy in the wavelength range from 2.0 to 4.7 μm[START_REF] Mumma | Organic Composition of C/1999 S4 (LINEAR): A Comet Formed Near Jupiter?[END_REF] have concluded that this comet has Jupiter. Nevertheless it is an Oort cloud comet, i.e. was expelled into the Oort cloud after formation. Comet C/1996 B2 (Hyakutake) has 1/a original = +0.002376 (MPEC 2001-X70) and probably has been in the inner solar system before. Comet WM 1 is in between the two others with 1/a original = +0.000533(MPEC 2000-N15).

(

  2001) observed a strong peak of the HCN production rate.A change of the particle size distribution leading to a strong enhancement of the number of small dust particles occurred inComet C/1996 B2 in mid-April 1996[START_REF] Fulle | Dust environment of Comet Hyakutake 1996B2[END_REF] constructed a dynamical model of this comet. They derived a change of the power index of the particle size distribution in this comet from ≈ -3.3 to < -4 in the interval 20-10 days before perihelion, i.e. before May 1.4 1996. The change was attributed to particle fragmentation, although other explanations may be possible. It occurred at a heliocentric distance of ≈ 0.55 AU. The formation of small dust grains was confirmed by the appearance of a silicate peak and an overheated dust coma in the thermal infrared spectral region[START_REF] Marsden | The Temporal Development of the Pre-perihelion Infrared Spectral A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT Energy Distribution of Comet Hyakutake (C/1996 B2)[END_REF]. This event was associated with

  1 . The data are plotted versus time. On top of the figure the corresponding heliocentric distances are shown. The brightness data of International Comet Quarterly show a small hump around Nov 19. The hump is real and not caused by a systematic effect of the bright moon affecting the visual magnitude estimates (Kammerer 2007, private communication).

  At the same time the method ensures positiveness of the flux

  All the different model solutions of Dec 03, despite of their significant differences, clearly indicate a gradual change of the particle size distribution toward a greater abundance of smaller dust grains as time progresses. This change is already suggested by the syndyne-synchrone plot of Dec 03 (Figure 3, center panel). It is possible, however, that our inversion method underestimates the contribution of large, recently emitted dust grains. Therefore the change to smaller particle sizes may not be as large as suggested by the models. (4) Total dust mass production (mass loss) ranges from 10 3 -10 6 g s -1 (Figure 14, center row). The dependence of dust mass loss on time differs for emission time orders n = 1 and n = 2, and, more important, the values derived from the three different observations are not consistent. The differences are caused by the method creating the inverse solutions. Most puzzling is the decrease of the mass loss with time (decreasing heliocentric distance) displayed by the model calculations of Nov 10 and particularly strong by those of Dec 03 (A similar case is decribed in inverse models of comet 46P/Wirtanen calculated by

Fig. 5 .

 5 Fig. 5. Terminal velocity (left scale) and ratio of radiative to gravitative force β (right scale) plotted versus particle size. Open squares and filled triangles denote

Fig. 10 .

 10 Fig. 10. The same as Fig. 9 but for Nov 19. The innermost level of intensity is 1.0 • 10 -6 in dimensionless pΦf units.

Fig. 11 .Fig. 13 .

 1113 Fig. 11. The same as Fig. 9 but for Dec 03. The innermost level of intensity is 3.2 • 10 -7 in dimensionless pΦf units.

  of the Max-Planck Institute for Aeronomy (now Max-Planck Institute for Solar System Research) was attached to the telescope. The focal reducer decreases the focal length of the 2m-telescope by a factor of 2.86, and at the same time increases the focal ratio of the telescope from f/8 to f/2.8.

Continuum images have been obtained with narrow-band filters centered at wavelengths 642 nm, 713 nm, and 853 nm. A CCD camera cooled with liquid nitrogen was used. The chip of this detector comprises 512 2 square pixels with

Table 1

 1 Observational circumstances. Exposure start-time. All exposures are 300 sec long. b heliocentric distance c geocentric distance d phase angle e angle between the comet-Earth vector and the comet orbit plane. Negative/positive angles refer to location of the Earth south/north of the orbital plane of the comet. f position angle of the prolonged Sun -Comet radius vector g Central wavelength and FWHM of filter used

	Date a	r b	Δ c	Scale	φ d	η e	α f	(λ/Δλ) g
	UT	AU	AU	km px -1 degree degree degree	nm/nm
	2001 Nov 10.79 1.553 0.615	443.5	18.5	-14.7	205.7 853.0/33.0
	2001 Nov 19.02 1.423 0.450	324.6	12.5	-2.5	157.0	641.6/2.6
	2001 Dec 03.70 1.186 0.317	228.6	44.6	+44.4	60.1	712.8/3.8

a

Table 2

 2 Distribution of the time bins in the Monte Carlo calculation (quasisynchrones) # True anomaly Day before perihelion Heliocentric distance r -1.5

	h	v 1	v 2

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

angular distribution is a sum of two terms: one term decribes particles emitted isotropically in all directions and the other an emission into a cone with half opening angle 30 • and with axis pointing to the Sun. But in our case, as described by eq. ( 7), for both terms the velocity of the particles smoothly depends on the angle between the emission position and the direction to the Sun. Thus, even in the case of isotropic emission (in the sense of an equal number of particles emitted in all directions) the velocity of the dust is reduced on the night side of the nucleus and drops to zero at the point opposite to the Sun. But in view of the inaccuracy of the model fit all the differences in the description of particle velocity with respect to Fulle are very likely insignificant. There is one important difference, however: In eq. ( 8) a 0 = 5 • 10 -5 cm and κ = 0.5 have been kept fixed all the time. Therefore we are unable to recognize fragmented small particles by their slow speed as compared to those having been accelerated by the outflowing gas (see [START_REF] Combi | The fragmentation of dust in the innermost comae of comets: Possible evidence from ground-based images[END_REF].

We use in this paper a different method for solving the inverse problem. After the detailed discussion of the achieved results, in particular after comparison of results belonging to different order combinations, a general assessment of the applied new inversion method seems appropriate. At first sight it may appear that the use of an inversion procedure as compared to trial-and-error is a mere convenience. Even now models using trial-and-error are proposed (Foster and Green, 2007). But the problem to determine dust size distributions and the temporal variation of the dust production rate from the shells M is nonunique, because of the finite ejection speed of the dust grains. Trial and error leaves the selection of the physically acceptable solution out of the manifold of solutions with similar χ 2 to chance. Nevertheless trial-and-error may be the method of choice if the data have a low signal/noise ratio (Lisse

from the nucleus or were created by fragmentation, but the latter appears much more likely. We would then favor a process like the one described by [START_REF] Boehnhardt | The possible role of organic polymers in the structure and fragmentation of dust in the coma of Comet P/Halley[END_REF] or [START_REF] Oberc | Small-scale dust structures in Halley's coma II. Disintegration of large dust bodies[END_REF], i.e. slow fragmentation caused by loss of an organic, polymeric glue.

In order to compare a dust mass production rate obtained from one of our models with the observed HCN and CN production rates we have added to [START_REF] Lara | Gas and dust in Comet C/2000 WM1 during its closest approach to Earth: Optical imaging and long-slit spectroscopy[END_REF] we get an approximate dust/gas ratio of 0.5. But in view of our discussion of section 5.4 this value may be erroneous by a wide margin.

Summary

In this paper we have analysed observations of the dust coma and tail of comet C/LINEAR (2000 WM 1 ) obtained on Nov 10, Nov 19 (appearance on an anti-tail) and Dec 03, 2001 with an inverse dynamical model. Particles with diameters d in the range 4.786 ≤ d ≤ 524.8μm, corresponding to β values 0.16 ≥ β ≥ 0.0016 were considered. As usual, the emission velocity of the dust particles from the nucleus was derived by trial and error. To produce inverse solutions (eq. 1) we have employed the method provided by [START_REF] Lemen | Coronal activity in F-, G-, and K-type stars. III -The coronal differential emission measure distribution of Capella, Sigma-squared CrB, and Procyon[END_REF] in a form generalized to two dimensions (emission time and particle size). This method enforces stability of the inverse process by imposing a varying degree of smoothness in the dependence of the flux vector F on emission time and

This certainly is an advantage. On the other hand the description of the flux vector F by Chebyshev polynomials (eqs. 14 and 15) adds another instability to the inversion process, when the method is used at higher orders. It is likely that in this case Lemen's method, in particular the two-dimensional version of this paper, preferentially selects solutions dominated by a low number of shells which then strongly contribute to the solution. To avoid this effect the method should not be used with orders higher than 3 in emission time or higher than 2 in particle size. Inversion methods based on maximum entropy considerations are perhaps preferable but they may not be as transparent and as easy to apply. In this paper we have for the first time presented different solutions of the inverse problem that describe the data to a very high degree equally well (Table 3,. For the first time the comparison of these solutions sheds light on the considerable ambiguity present in the process of inversion

of an improperly posed problem. The dynamical modeling of cometary dust tails may perhaps be the best way to derive particle properties from remote sensing observations (Fulle, 2004), but even the best way still leaves much to be desired. All other inversion methods available in the literature and mentioned in Section 5.5 also contain parameters that allow to produce solutions differing in the amount of smoothing or regularization. It is very important that other researchers attempt to analyze their inversion methods and investigate the selection process of the "well behaved" fitting solution with respect to the set of available solutions of the non-unique, ill-posed, problem.

Acknowledgments

Irina Kulyk provided valuable assistance at the telescope. We are grateful to Svetozar Zhekov for drawing our attention to Lemen's method [START_REF] Lemen | Coronal activity in F-, G-, and K-type stars. III -The coronal differential emission measure distribution of Capella, Sigma-squared CrB, and Procyon[END_REF] and for discussions concerning this subject. We thank Nikolaj Kiselev and Viktor Tishkovets for fruitful discussions. Brian Marsden communicated the values of the original inverse semimajor axes of the comets mentioned in Section 5.6. We thank the referees J. 9)) is plotted versus quasisynchrone number (see Fig. 4).

In the first column (labeled "Sum") on the right side of the velocity plot the sum of all shells further to the right, i.e. an image of the corresponding quasisynchrone, is

shown. The rest of the columns presents images of 11 of the 51 shells that form the quasisynchrone. Their β-values are indicated in the bottom. In the lower rows some fields belonging to more strongly accelerated shells are empty because the shells have left the field of view. When the ejection velocity u increases, the shells and quasisynchrones become wider. This effect is reduced in the most recently emitted shells as in those there was not enough time for expansion. See Table 2 for the description of quasisynchrone parameters. [START_REF] Biver | Radio wavelength molecular observations of comets C/1999 T1 (McNaught-Hartley), C/2001 A2 (LINEAR), C/2000 WM 1 (LINEAR) and 153[END_REF][START_REF] Lara | Gas and dust in Comet C/2000 WM1 during its closest approach to Earth: Optical imaging and long-slit spectroscopy[END_REF][START_REF] Schleicher | Photometry and Imaging of Comet C/2000 Wm1 (Linear)[END_REF] and our maximum and minimum values of the dust mass production rate from Nov 19, order combination (2,1). The horizontal parts of the steps indicate the time interval of a single quasisynchrone. 5 mag (right ordinate) correspond to a factor 100 in the production rate (left ordinate). Therefore the magnitude and production rate slopes are directly comparable. See text.