text:
# Realistic error estimates on kinematic parameters

## 1 Introduction

The bulk of kinematic information available on elliptical galaxies is nowadays retrieved through observation and analysis of the line-of-sight velocity distributions (LOSVD’s) of the stars in these galaxies.

The basic idea behind the study of LOSVD’s is that a galaxy spectrum can be generated by the spectrum of a star of similar type that has been redshifted by the rotation of the galaxy and smeared out by the velocity dispersion of the stars (Minkowski, 1962).

The study of LOSVD’s as we know it now started about 30 years ago, when near the end of the seventies, Sargent et al, (1977) and Tonry & Davis, (1979) published surveys of galaxy redshifts and velocity distributions. The methods they presented in their papers were widely used afterward:

* Fourier quotient method (Sargent et al,, 1977): The LOSVD is obtained by dividing the Fourier transform of the galaxy spectrum by the Fourier transform of the template spectrum, and transform the result back into pixel space. Although it seems to be natural to use Fourier space, this has also some disadvantages. The error calculation becomes complicated due to the correlations in the quotient and the method is very sensitive to template mismatch.
* Cross-correlation method (Tonry & Davis,, 1979): The peak in the galaxy-template correlation function is fit by a Gaussian function. The method also provides a way to obtain error estimates on the kinematic parameters.

These techniques assume a Gaussian form for the broadening functions for mainly practical reasons: it provides a way of filtering the noise in the spectra and it allows to calculate convolutions analytically. The mean and variance determining the Gaussian function were considered as parameters describing the mean streaming and stellar velocity dispersion.

About ten to fifteen years later, the quality of the data had improved a lot and new techniques were developed trying to get more information out of the observed spectra than just two parameters. The main improvement was the elimination of the assumption of a Gaussian form for the LOSVD. There are mainly two criteria that can be used to classify these techniques: (a) whether a non-parametric LOSVD or a parameterized LOSVD is derived. In case of high S/N data, a non-parametric LOSVD yields all information you can get out of a spectrum. But because of the presence of noise there is always a need for some sort of smoothing. There are several ways to do this, either by adapting some filter or by using a smoothing parameter. For these non-parametric fitting methods there is not always a simple way to come to an error estimate. (b) whether the fit is performed in Fourier space or in pixel space. A treatment of the problem in pixel space has the advantage that parts of the spectrum can be easily eliminated and allows an easier error estimate.

A lot of methods have been proposed:

* Fourier fitting method (Franx et al.,, 1989): (parametric in Fourier space) A convolution of the LOSVD with the template is fitted to the galaxy spectrum in Fourier space. Error estimates are derived from a least squares fit.
* Fourier correlation quotient method (Bender,, 1990): (non-parametric in Fourier space) The method is based on the deconvolution of the peak of the template-galaxy correlation function with the peak of the autocorrelation of the template star. This approach is shown to be less sensitive to template mismatch than the standard Fourier quotient method. The deconvolution itself is done in Fourier space. Noise filtering is done by applying a Wiener filter. Error estimates can be retrieved from Monte Carlo simulations.
* Direct fitting method (Rix & White,, 1992): (non-parametric in pixel space) The galaxy spectrum is assumed to be a superposition of various components: (1) continuum terms, (2) a set of broadened stellar spectra where the broadening function and the stellar spectra can vary, (3) a noise component. The best combination of components and the errors on the results are obtained from a least squares minimization.
* unresolved Gaussian deconvolution (Kuijken & Merrifield,, 1993): (non-parametric in pixel space) The LOSVD is modeled as a sum of a set of Gaussian distributions that are uniformly spaced in velocity. The best fit is determined by means of a least square minimization, which also serves as basis for error determination.
* Bayesian method (Saha & Williams,, 1994): (non-parametric or parametric in pixel space) The concept is to take a trial LOSVD and convolve it with a template spectrum into a model spectrum. If this model is the real one, the difference between the model spectrum and the observed spectrum must be the noise. With probability theory it can be calculated what the probability is that the trial LOSVD is the true one, given the data. Confidence intervals for broadening functions are obtained from Monte Carlo simulations.
* Cross-correlation method update (Statler,, 1995): (parametric in pixel space) This is a generalization of the galaxy-template cross-correlation method to arbitrary parameterized line profiles, where the correlation function is fitted by the model line profile using a least squares technique. For a well resolved galaxy spectrum, the result is fairly insensitive to template mismatch. Errors follow out of standard statistics involved in the fitting technique.
* penalized likelihood (Merritt,, 1997): (non-parametric in pixel space) The method is based on a penalty function that is large for any noisy function and zero for a Gaussian function. This method allows to fit high quality data very accurately, whereas also very noisy data can be fitted, in the limit by a simple Gaussian. Monte Carlo simulations are performed to come to an error estimate.

When parameters based on the integrated moments of the LOSVD are sought for, one can also obtain these parameters in a one step process. In this case, a convolution of a parameterized expression for the LOSVD with a stellar template spectrum is matched to the observed galaxy spectrum.

* Gauss-Hermite fitting method (van der Marel & Franx,, 1993): (parametric in Fourier space) A Gauss-Hermite parameterization is used for the LOSVD, a least squares minimization is performed in Fourier space.
* Gauss-Hermite direct fitting method (van der Marel et al.,, 1994): (parametric in pixel space) The LOSVD is written as a truncated Gauss-Hermite series and the least squares fit is performed in pixel space.

It is generally believed that, if the data are good enough to apply non-parametric fitting, such an approach implies less bias than a parametric fit. However, having a non-parametric LOSVD, it is still very useful to describe the obtained profile by just a few characteristic numbers. Higher order information can be obtained either by calculating the moments of the velocity distribution or by fitting a parameterized function to it. Since there is a real danger of confusion when dealing with the moments of the LOSVD, we define them here explicitly:

$$v_p=𝑑v_pv_p(v_p),$$
(1)
and

$$\mu _k=𝑑v_p(v_pv_p)^k(v_p),k1,$$
(2)
with $`(v_p)`$ the LOSVD. The kinematic moments are the mean projected rotation $`v_p`$, the projected velocity dispersion $`\sigma _p^2=\mu _2`$, the skewness $`\zeta _3=\mu _3/\sigma _p^3`$ and the kurtosis $`\zeta _4=\mu _4/\sigma _p^4`$.

The calculation of higher order moments out of LOSVD’s may be troublesome because of the strong dependence of the moments on the wings of the profiles. These wings unfortunately are the hardest to retrieve from observations. Therefore, in most practical cases, a parameterized function is fitted to the LOSVD. A truncated Gauss-Hermite series is often used (Gerhard, 1993,van der Marel & Franx, 1993):

$$(v_p)=\gamma \mathrm{exp}\left[\frac{1}{2}\frac{(v_pv)^2}{\sigma ^2}\right]\times \left[1+\underset{j=3}{\overset{n}{}}h_jH_j\left(\frac{v_pv}{\sigma }\right)\right].$$
(3)
In this expression, $`v`$ and $`\sigma `$ are free parameters, they are also called ’mean velocity’ and ’velocity dispersion’. The functions $`H_j(x)`$ are the $`j`$-th order Gauss-Hermite polynomials, and their coefficients $`h_j`$ are regarded as the higher order parameters. The coefficient $`h_3`$ is the third order shape parameter and is an indication for anti-symmetric deviations from a pure Gaussian function in the form of leading or trailing tails. The coefficient $`h_4`$ is the fourth order shape parameter and is an indication for symmetric deviations from a pure Gaussian LOSVD in the form of broader or smaller wings.

Kinematic information up (to parameters corresponding) to the fourth order moment of the LOSVD is becoming available for a large sample of ellipticals, see e.g Bender et al. (1994), Kronawitter et al., (2000), Halliday et al., (2001). These data sets are a prerequisite for the construction of dynamical models of elliptical galaxies. In some cases it can be noticed that the error estimates on the parameters are smaller than the scatter among neighbouring data points in the kinematic profiles. This is especially striking for the so-called higher order moments (third and fourth order moment). Examples of this behaviour can be seen in some of the data presented by Bender et al. (1994), Kronawitter et al., (2000), Halliday et al., (2001). A possible interpretation of this fact is that the real errors on the kinematic parameters are larger than the error estimates given by the parameter fitting algorithms.

The derivation of these error estimates relies on assumptions that are generally not completely met. Though people are aware of this, it is generally unclear how large the discrepancy between reality and assumptions is.

One can come up with at least two good reasons why it is important to have realistic error estimates for dynamical modelling. First, for the construction of dynamical models, it is convenient to work with a goodness of fit indicator, which is meant to be an objective indicator of how good the model fits the data. Such goodness of fit indicators need error estimates on the data. However, if these are not realistic errors, it is not straightforward to interpret the values of this goodness of fit indicator, let alone derive meaningful confidence intervals for the models. In some cases, awareness of this problem may prompt one to adopt modified error estimates for use in the goodness of fit indicator (e.g. Kochanek, (1994), Cretton et al. (2000), Gerhard et al., (1998), Saglia et al., (2000), Kronawitter et al., (2000)). However, playing with the error estimates on individual data points is equivalent to using different relative weights for the data points, and this is likely to complicate the interpretation of a goodness of fit indicator rather than simplify it.

A second issue of concern is that these parameters are often used to determine a realistic dynamical model and to constrain the potential of a galaxy at the same time. In particular the velocity dispersion and the fourth order parameter play an important role in inferring the existence of central black holes or dark matter haloes around elliptical galaxies. In practice this is done by assuming a potential for the system and trying to find a dynamical model based upon the selected potential that simultaneously fits the velocity dispersion and fourth order parameter within the given error estimates. If no such model can be found, it is concluded that the assumed potential is not the right one, and the amount of dark matter is changed. If the error estimates are smaller than the real errors and the conclusions on the potential critically depend on them (e.g. Rix et al., (1997), Kronawitter et al., (2000), Cretton & van den Bosch, (1999)), some potential-model combinations may be rejected wrongly. As an inevitable consequence, it is possible that in some cases the use of realistic error estimates in dynamical modelling will weaken the dynamical evidence for dark matter in elliptical galaxies.

The method to derive realistic error estimates on kinematic parameters that is presented in this paper reaches out toward a need on the observational side and the modelling side. The method makes an implicit diagnosis of the noise distribution and uses the characteristics of this distribution as the basis of a new error estimate.

In section 2, a diagnosis of the problem with error estimates is given. An outline of the new method is presented in section 3, this is elaborated and illustrated in section 4. A discussion and conclusions can be found in sections 5 and 6 respectively.

## 2 Diagnosis of the problem

In the absence of external sources of errors, a raw spectral image contains data points that follow a Poisson distribution, the variance of which is denoted by ”Poisson noise”. Another characteristic of a raw spectral image is that the data points in the image are independent. Hence, the power spectrum of the raw signal would match the definition of white noise: it has constant power in the frequency domain.

Before kinematic parameters can be calculated, the spectra have to be brought in a state that is referred to as ’cleaned’ and ’calibrated’. Bringing the spectra in such a condition requires a number of image processing steps, the so-called ’standard reduction techniques’, for which adequate tools are available in image analysis packages such as MIDAS or IRAF. But besides bringing the observed spectra in good condition, these reduction steps also have a non-negligible influence on the noise distribution.

For example, the removal of bad pixels or cosmic ray events, where the values of two or more adjacent pixels are replaced by an average over a surrounding region, clearly introduce dependencies in the data. Likewise, dependencies are introduced by any operation where some kind of interpolation is involved, like e.g. calibration. Moreover, a sky subtraction does certainly not guarantee a perfect elimination of the sky in the data, regardless the level of sophistication that is used to perform the operation. The residual of this correction also adds up to the noise distribution. In some cases the above effects can be the source of considerable uncertainties on the derived kinematic parameters, while this is not reflected by error estimates based on Poisson noise.

In many cases, a least squares technique is used to derive the kinematic parameters and error estimates from the data. The statistical interpretation of this method relies on the assumption that (1) the noise is independent and (2) Gaussian distributed on the input data. For large numbers, a Gaussian distribution is a fairly good approximation for a Poisson distribution and therefore the second assumption is valid. It is clear from the above arguments that the first condition generally is not met after data reduction. As a consequence, the errors derived from standard statistics will in many cases differ from the real errors on the kinematic parameters.

Sometimes, Monte Carlo simulations are used (Bender et al. (1994), Statler, (1995)) to estimate the uncertainties on the derived parameters. For the realization of synthetic galaxy spectra, a Gaussian noise distribution with given S/N is used. So also this method may show a similar tendency to underestimate the errors.

Accepting the fact that the cleaned and calibrated spectral images come as they are, it is possible to obtain error estimates that are more realistic estimates than those currently found in literature. This is the purpose of the present paper.

## 3 Outline of the method

To get an idea of more realistic errors on the parameters and to overcome the above problems we propose the following scheme. The steps are described in more detail in the following sections.

* Perform a fit of a model galaxy spectrum to the observed galaxy spectrum and determine the residual of this fit. The model galaxy spectrum is a convolution of the selected template spectrum with a LOSVD. For this step it is important to use an expression for the LOSVD that offers enough freedom. The purpose is to fit every part of the galaxy spectrum that can reasonably be fit by a convolution of the template spectrum with a general and smooth LOSVD.
* The residual of this fit contains (1) the Gaussian noise on the raw spectrum, (2) the features left or introduced by several reduction steps. The errors estimated with the proposed method take into account artifacts from flatfielding, cosmic remnants, effects of rebinning and remnants of sky removal, for as far as the latter has been removed reasonably well. On the other hand, the method does not account for errors in continuum subtraction nor template mismatch. Hence, the residual of the fit can be considered as a realization of the ”observed noise” involved in the problem: the term ”observed noise” applies to ’cleaned and calibrated’ spectra and is used to make a distinction with Poisson noise. Where pure white noise is featureless and transforms into a flat power spectrum, the power spectrum of the residual will be different. A smooth representation of the power spectrum of the residual will be taken as sufficient information on the main characteristics of the observed noise.
* The next step is to generate a number of synthetic galaxy spectra, that have a noise distribution with the same power spectrum as the observed noise, and to use these synthetic spectra to determine the spread on the kinematic parameters. These equivalent noise distributions that carry the characteristics of the original observed noise can be obtained by multiplying the Fourier transform of a white noise profile with a smooth representation of the power spectrum of the residual, and transforming this back.

The idea of using Monte Carlo simulations to determine error estimates of course is not new, but the main point is that in the present method one tries to model the real characteristics of the noise.

A demonstration of the method on realistic galaxy spectra is given in section 4.

## 4 NGC 3258

### 4.1 Observations

NGC 3258 is an E1 galaxy and was observed with the ESO-NTT telescope in the nights of 27-28/2/2001 (64.N-192). Spectra of the major axis were taken using the red arm of EMMI, covering the Ca II triplet. Grating #7 was used, having a dispersion of 0.66 Å/pix. The detector was a Tektronix CCD with 2048$`\times `$2047 pixels, 24$`\mu m\times `$24$`\mu m`$ in size and with a pixel scale of $`0.27^{\prime \prime }`$/pixel. A slit width of $`1.5^{\prime \prime }`$ thus yields a spectral resolution of 3.67 Å FWHM, resulting in an instrumental dispersion of about 54 km/s in the region of the Ca II triplet. For NGC 3258, several exposures of 3600 sec were taken (in total 11 hours). A number of standard stars (G dwarfs and K and M giants) were also observed.

#### 4.1.1 Standard data reduction steps

Standard reduction steps were applied to these spectra. Out of the bias and dark exposures a correction term was derived. Several domeflat images were used to come to an appropriate flatfield image. Cosmics were removed using a top hat filter. After filtering, the images were inspected carefully and remaining cosmics were also removed by hand.

For the wavelength calibration, lamp spectra were taken just before or after each of the spectroscopic observations of the galaxy or the stellar template stars. For each row in the lamp spectra images, a polynomial was constructed to transform the pixel scale into wavelength scale. The same coordinate transformation was then applied to the galaxy or star spectra. The calibrated images have a step of $`0.3`$Å.

For the airmass correction, the mean value of the airmass at the beginning and end of the exposure was used. The contribution of the sky to the spectra was estimated from an upper and lower region of the image, where there was no contribution of the light of the galaxy or template star. Several galaxy spectra were reduced separately, aligned and combined into one galaxy spectra image.

Spatial rebinning of the galaxy spectra resulted in data with a S/N $``$ 50/bin (about 70 in the centre). Only a small part of the spectrum, from 8500 Å to 8750 Å was used to determine the kinematic parameters.

### 4.2 The observed noise on the data

The idea is that the residual of the observed galaxy spectrum with a sufficiently flexible model spectrum can be attributed to some sort of noise involved in the problem. For a given stellar template spectrum it is important to consider a representation for the LOSVD that offers enough degrees of freedom to be able to find the best fitting model spectrum.

For this reason, nonparametric LOSVD’s in combination with a flexible method to generate them are chosen to perform this fit. We have also chosen not to restrict ourselves to positive LOSVDs, because this makes the fit more general and at this point there are no physical restrictions on the LOSVDs. The method adopted here performs a least squares fit, using a set of cubic spline basis functions to represent the LOSVD. An outline of the method can be found in Appendix A. However, for this step, any non-parametric fitting method with enough degrees of freedom will do. The model spectrum was a linear combination of convolutions with a number of observed template stars in which the coefficients were chosen on the basis of a goodness of fit estimator.

The non-parametric fit is only used to get an idea of the fraction of the signal in the galaxy spectrum that cannot be accounted for by a model spectrum with the chosen template mix, and hence that will be part of the residual. Furthermore, also traces of sky lines, or incompletely removed cosmics will contribute to the residual of this fit.

As an illustration, the results of the method for two different positions along the major axis of NGC 3258 are shown in figure 1 for $`r=0.8^{\prime \prime }`$ and figure 2 for $`r=11^{\prime \prime }`$. The upper panel displays the LOSVD, obtained from a non-parametric fit with 10 degrees of freedom. The middle panel shows the observed galaxy spectrum and the fit. The lower panel shows the residual of the fit. From these lower panels is it clear that the residual at $`11^{\prime \prime }`$ (values mostly between -100 and 100) is much larger than the residual at $`0.8^{\prime \prime }`$ (values mostly between -60 and 60). This immediately shows that the Poisson noise is not the only source of errors. The galaxy spectrum at $`0.8^{\prime \prime }`$ has much more counts in absolute value than the galaxy spectrum at $`11^{\prime \prime }`$, hence pure Poisson noise is higher at $`0.8^{\prime \prime }`$ than at $`11^{\prime \prime }`$. The surplus non-Poissonian contribution to the residual at $`11^{\prime \prime }`$ is likely to come mostly from an imperfect sky subtraction.

### 4.3 The power spectrum of the residual

To see what the magnitudes of the various frequency components in these residuals are, their power spectra have been calculated. They are presented in figure 3, the residual at $`0.8^{\prime \prime }`$ is in the upper panel, the residual at $`11^{\prime \prime }`$ in the lower panel. The scale on the horizontal axis is the period of the power spectrum, instead of the more widely used frequency. This period (labelled $`\mathrm{\Delta }\lambda `$) is directly related to the wavelength scales of the features in the original spectrum. Because of the wide range in values for $`\mathrm{\Delta }\lambda `$, a logarithmic axis is plotted. On top of the residual power spectra in figure 3, a sample power spectrum for the expected white noise at these radii is presented. There is a clear difference between both profiles.

Having these power spectra representations, one would like to know which part is playing an important role for the determination of real errors on the kinematics. In other words, which frequencies can have an important impact on the derived kinematic parameters?

This can be investigated by means of some simulations. Galaxy spectra with various noise characteristics can be created starting from a galaxy spectrum obtained after convolution of a template spectrum with a known LOSVD and adding wave functions (sine and cosine functions) with different frequencies. For these artificial spectra, a best fitting model spectrum was determined using a least-squares minimization. The model spectrum was composed with the template spectrum and a LOSVD expressed as a truncated Gauss-Hermite series. This analysis of the artificial spectra gives new values for the kinematic parameters, that will differ from the original kinematic parameters. The results of such simulations as a function of the scale size of the waves can be seen in figure 4 (upper panels for d$`v`$/$`\sigma `$ and d$`\sigma `$/$`\sigma `$, lower panels for d$`h_3`$ and d$`h_4`$).

It is clear that the measurement of these kinematic parameters is sensitive to superpositions of wave functions with wavelengths ($`\mathrm{\Delta }\lambda `$) between 10 Å and 400 Å (hashed regions in figure 4). This $`\mathrm{\Delta }\lambda `$ range is independent of the wavelength range used in the analysis, but scales with $`\sigma `$ (which is 325 km/s in this case).

This learns that the region of interest in the power spectrum of the residuals lies between $`\mathrm{\Delta }\lambda =10`$Å and $`\mathrm{\Delta }\lambda =400`$Å. It is also there (see figure 3) that the residual power spectra overshoot the white noise.

For short wavelengths, the residual power spectrum seems to have lower values than the Poisson noise power spectrum. At first sight this seems to be puzzling. However, this is a result of the interpolations that took place in the course of the data reduction process. Interpolation always implies some sort of smoothing, which means that features with a small wavelength scale disappear from the image.

Moreover, oversampling lies at the origin of the transfer of an amount of noise coming from small wavelength scale to larger wavelength scale. The stronger the oversampling, which has the artificially enlarging of feature wavelength scales as a result, the stronger this migration effect. In the case of NGC 3258, the original spectrum had a step of $`0.66`$Å/pixel, while the spectra after calibration have a step of $`0.3`$Å/pixel.

### 4.4 Realistic error estimates

The power spectra of the residuals clearly show that the assumption of independent (white) noise is not a valid one. This implies that the classical statistical tools for calculating error estimates out of least squares fitting techniques cannot be used.

The only option that is left, is to determine the uncertainties on the parameters through Monte Carlo simulations. In this case it is important to work with simulated galaxy spectra that have the same noise characteristics as the original galaxy spectrum. Our method proposes to achieve this by using simulated galaxy spectra that give residuals that follow the same power spectrum as the initially determined residual.

To realize this, a smooth representation of the residual power spectrum is created. This is illustrated in figure 3, by means of the dashed lines. In this case, each pixel was replaced by a weighted mean over 20 pixels at each side taking a linearly declining weight function into account. This smooth function has basically the same behaviour as the power spectrum, and hence carries the information that distinguishes the observed noise distribution from white noise. New realizations of this noise distribution can be calculated by multiplying a white noise spectrum with this fitted function. In this way, the simulated galaxy spectra for the Monte Carlo simulation are constructed.

To these spectra, a model spectrum composed with the template spectrum and a parameterized expression for the LOSVD is fitted. For the parameterization, a Gauss-Hermite decomposition was used, following van der Marel & Franx, (1993). The best fitting values for $`v`$, $`\sigma `$, $`h_3`$ and $`h_4`$ were obtained using the Levenberg-Marquardt method as described in Press et al., (1989).

The error estimates are the root mean squares of the kinematic parameters that were derived from 100 realizations of the galaxy spectra. The results are shown in figure 5, from upper to lower panel: the mean rotation velocity, the velocity dispersion, $`h_4`$ and $`h_3`$. The vertical error estimates present realistic errors. The horizontal lines correspond to error estimates based on least squares statistics.

The mean velocity reaches a maximum value of about 56 km/s and the realistic error estimates lie between 5.7 km/s and 11.3 km/s (with a mean error of 7.5 km/s and a mean error from least squares statistics of 4 km/s). The central velocity dispersion is 344 km/s, with an error of 8 km/s. The mean error on $`\sigma `$ is 10.4 km/s, the mean error on $`\sigma `$ based on least squares statistics is 3.8 km/s. For large radii, the profile is antisymmetric.

The profile for $`h_3`$ goes through the origin, apparently there is no serious problem with template mismatch. The mean error on $`h_3`$ is 0.02, the lowest value is about 0.013, the highest value is 0.039 (measured in an outer point). From least squares statistics, the mean error on $`h_3`$ is 0.01.

The $`h_4`$ profile is approximately symmetric, with low values in the centre and higher values for the outer data points. The central data points show quite some scatter, but this is not inconsistent with the errors (at least 0.014) on the data in that region. The errors increase outwards, the outer data points have an error of 0.049 (positive side) and 0.03 (negative side). The mean error on $`h_4`$ is 0.023, the mean error on $`h_4`$ following least squares statistics is 0.01.

For the $`h_3`$ and $`h_4`$ profiles, we compared the data points at either side of the centre. The amount of scatter in the profiles can be expressed as the mean distance between corresponding data points at positive and negative radius, $`\sqrt{(_{i=1,n}(|y_i||y_i^{}|)^2)/n}`$, with $`n`$ the number of data points at positive or negative radius, and $`y_i`$ and $`y_i^{}`$ the values for the parameters in corresponding points at positive and negative side of the centre. Fore $`h_3`$ the value of this scatter indicator is 0.019 and for $`h_4`$ this is 0.031. This means that for both parameters, the mean errors derived with the proposed method are close to this scatter indicator, while the mean errors following from least squares statistics are clearly smaller.

## 5 Discussion

To verify and validate the algorithm and its implementation, we now explore a few situations in which (following the discussion in section 2) conventional and more realistic error bars based on the method described in this paper are expected to give different results by means of simulations. Using simulations has the advantage that also the true errors can be calculated and hence the error estimates can be compared.

### 5.1 Impulse noise

A fraction of the noise with non-Poissonian character is manifested as spikes that overshoot the Poisson noise and that are relatively small on wavelength scale. This situation can be mimicked in a simulation by adding impulse noise to a spectrum with otherwise pure Poisson noise.

Synthetic model spectra were obtained by convolving a LOSVD characterized by $`v=h_3=h_4=0`$ and $`\sigma =250`$ km/s with a K3 giant and adding Poisson noise. Different types of model spectra were then created by contaminating 10 $`\%`$ of the pixels in the spectrum with impulse noise with an amplitude of 3, 6, 9, 12 or 15 times the amplitude of the Poisson noise.

The RMS deviation on the kinematic parameters calculated out of 100 equivalent realization of one type of synthetic spectrum are considered to be the true errors on the kinematic parameters in that specific situation. For each type of spectrum, also the analyses following the scheme outlined in section 3 was applied (these results will be referred to as realistic errors). Likewise, conventional errors were calculated with a standard least squares technique. These results are presented in figure 6. The true errors are indicated in asterisks, the realistic errors in squares and the conventional errors in diamonds.

It is clear that the conventional errors returned by a standard technique are independent of the degree of contamination and can underestimate the true errors considerably. For an impulse noise with an amplitude of 3 times that of the Poisson noise the conventional error underestimates the true error by $`27\%`$ for $`v`$, $`37\%`$ for $`\sigma `$, $`31\%`$ for $`h_3`$ and $`35\%`$ for $`h_4`$. The realistic errors on the other hand follow the true errors much better. For the same situation the difference between realistic errors and true errors is $`10\%`$ for $`v`$, $`5\%`$ for $`\sigma `$, $`0\%`$ for $`h_3`$ and $`14\%`$ for $`h_4`$. For serious contamination, with an impulse noise with an amplitude of 15 times that of the Poisson noise, the realistic error bars seem to slightly overestimate the real errors, but this is hopefully no longer a realistic situation.

### 5.2 Sampling

Changing the original sampling step of a spectrum has an impact on the error propagation that is completely neglected when Poisson noise is used as basic assumption for the error calculation.

For the simulations, a number of synthetic spectra were obtained with the same LOSVD and template star as above, but the number of pixels was different for each spectrum. In the end, all synthetic spectra were resampled in order to end up with spectra with $`N`$ pixels. This means that for the final spectra, the wavelength scale was $`n`$ times the original wavelength scale, with $`n`$ equal to 1.5, 2, 3 or 4, in order to mimic oversampling, and also $`n`$ equal to 0.5 and 0.25, thus mimicking undersampling.

These sampling factors are chosen only for the sake of illustrating the method. These rebinned spectra were analysed using the technique described in section 3, yielding realistic errors, and using a standard technique for the conventional errors. Again, the real errors were calculated out of Monte Carlo simulations with 100 spectra. The resulting error estimates are presented in figure 7, in squares for the realistic errors and diamonds for the conventional errors. The true errors are indicated in asterisks. Again, the conventional errors are independent of any change in amount of over- or undersampling, while the realistic errors mainly follow the true errors. The change in errors scales roughly with $`\sqrt{n}`$, with $`n`$ the sampling factor.

This result shows that care should be taken if kinematic parameters and conventional errors are estimated from spectra that are logarithmically rebinned. In that case the wavelength scale is clearly changed, but in such a way that some regions are more compressed and other regions are stretched.

### 5.3 Template mismatch

We do not claim that the proposed method is able to cope with the errors coming from template mismatch. But also in this case, the errors estimated with the method presented in this paper are closer to the true errors than the conventional errors obtained by standard techniques. The synthetic spectrum, created with a K3 giant, was also analysed with the following templates: G2 dwarf, G5 dwarf, K1 giant, K4 giant, M1 giant, M2 giant. The results are presented in figure 8. The horizontal axis shows the sequence of template stars, the bars indicate the errors: black for the true error, dark gray for the conventional error and light gray for the realistic error (indicating that the error is obtained using the technique presented in this paper). The results for $`\sigma `$ are in the upper panel, the results for $`h_3`$ in the middle panel and the results for $`h_4`$ in the lower panel. The results for $`v`$ are not shown because the correction for the radial velocity of the template star would introduce an additional uncertainty.

It is clear that the true errors (black bars) are larger than the conventional errors and the realistic errors. For all three kinematic parameters, it is remarkable that the conventional errors show little variation if different stellar templates are used, whereas the differences in values for the ”realistic” errors are larger. Moreover, in cases with large true errors, the ”realistic” errors are closer to these values than the conventional errors. It seems that the errors obtained with this method can be used as indicators for template mismatch in the sense that larger errors indicate a poorer fitting template, whereas the conventional errors clearly cannot be used as such.

## 6 Conclusions

As the technical limitations on the quality of observed data are diminishing, the numerical signal processing starts putting limitations on the information content of the observations. Therefore, it is important to have realistic estimates of the errors on the data.

In this paper, we propose a new method to calculate error estimates on kinematic parameters derived from spectroscopic observations. The starting point is the realization that the data points in the ’reduced and calibrated’ spectral images do not follow a Poisson distribution and are not independent from pixel to pixel. This implies that the error estimates derived with classical tools (Monte Carlo simulations with Gaussian noise or statistical interpretation of least squares techniques) will underestimate the true errors in many cases.

The proposed method first determines the characteristics of the real noise on the data. For this, a flexible fitting method offering enough degrees of freedom is used. The residual of the fit is considered as observed noise. It is shown that the power spectrum of this observed noise can be very different from white noise. Realistic error estimates are obtained through Monte Carlo simulations with synthetic galaxy spectra, where the noise distribution follows the same power spectrum as the observed noise. In this way, errors are estimated taking the real characteristics of the noise into account, instead of relying on a Gaussian (or Poisson) noise distribution. Moreover, the proposed scheme can be applied in combination with several methods that are currently used to derive kinematic parameters from spectroscopic data.

The method was tested on spectroscopic observations of NGC 3258. In this particular case, the realistic errors are almost a factor of 2 larger than the errors based on least squares statistics. What is most important, is that the realistic error estimates are more consistent with the scatter among neighbouring data points in the kinematic profiles.

Simulations with spectra containing Poisson noise and impulse noise confirm that the proposed method offers error estimates that are closer to the true errors than conventional error estimates. Moreover, from simulations with synthetic spectra it becomes clear that an oversampling of the spectra results in an underestimate of the errors on kinematic parameters when simple least squares statistics are used. This may be an important consideration when kinematic parameters are calculated from logarithmically rebinned spectra. Although this method is not able to calculate realistic errors on kinematic parameters obtained with an ill matching template star, the realistic errors do trace the situation of template mismatch, in the sense that they become larger with increasing template mismatch, unlike the usual error estimates.

## Acknowledgments

## Appendix A Cubic spline fit

In this paper, a new method is used to perform a non-parametric fit of a model spectrum to an observed galaxy spectrum. The method performs a least squares fit, using a set of cubic spline basis functions to represent the LOSVD.

Each basis function is a third-degree cubic spline. There are two types of cubic spline basis functions, each a composition of two polynomials, as illustrated in figure 9.

In this figure, the cubic splines are defined in the interval \[-1,1\]. The first third degree basis function (left in figure 9) is a composition of polynomials I and II. The second basis function (right in figure 9) is a composition of polynomials III and IV. These polynomials are:

$`I`$ $`=`$ $`13x^22x^3=(1+x)^2(12x)`$ (4)
$`II`$ $`=`$ $`13x^2+2x^3=(1x)^2(1+2x)`$ (5)
$`III`$ $`=`$ $`x+2x^2+x^3=(1+x)^2x`$ (6)
$`IV`$ $`=`$ $`x2x^2+x^3=(1x)^2x`$ (7)
Function I is defined in the interval  and meets the requirements that $`y(1)=0=y^{}(1)`$, $`y(0)=1`$, $`y^{}(0)=0`$, (prime denoting first order derivative). This behaviour is mirrored with respect to the vertical axis into the region \[-1,0\], resulting in function II. The composition of I and II gives a symmetric cubic spline basis function. The first two conditions assure that the function declines smoothly to 0 at the edges. The other two conditions give the profile a flat peak on the y-axis.

Function IV is defined in the interval  and meets the requirements that $`y(1)=0=y^{}(1)`$, $`y(0)=0`$, $`y^{}(0)=1`$. This behaviour is mirrored with respect to the centre into the region \[-1,0\], resulting in function III. The composition of III and IV is an antisymmetric cubic spline basis function. Again, the first two conditions assure that the function declines smoothly to 0 at the edges. The other conditions make sure the function goes through the centre and give a steeply declining or inclining profile near the centre.

For the practical implementation, the interval where the fit is performed is divided in $`(n1)`$ sub-intervals, using $`p_i`$, $`i=1,\mathrm{},n`$ points. This is illustrated in figure 10. Each triplet of adjacent points $`p_{i1},p_i,p_{i+1}`$ is used to define a set of cubic splines that is included in the fit.

The choice of the number of sub-intervals or basis functions is free. The less sub-intervals are considered, the more the line profile is smoothed. There is no a priori criterion that can be used to decide how many of these sub-intervals have to be used to obtain a good fit. Instead, a number of fits with an increasing number of sub-intervals (hence basis functions) is performed and the $`\chi ^2`$ values of the fits then indicate when the number of cubic spline pairs offers enough degrees of freedom to come to a good fit. This is illustrated in figure 11; once a sufficient degree of freedom is reached, the values for $`\chi ^2`$ start decreasing only very slowly when the number of nodes is increased. If $`n`$ points are used, $`(n1)`$ sub-intervals are used for the fit, yielding $`2n4`$ degrees of freedom.

Kuijken & Merrifield (1993) came up with a similar idea, but they used Gaussian functions for the decomposition of the LOSVD. The cubic splines used here are for the same degrees of freedom slightly more flexible in fitting. The fit is not restricted to strictly positive results. If one wants to obtain physical solutions with this method, linear constraints can be added to the $`\chi ^2`$ fit.


format: text

text:
# Multisymplectic Geometry Method for Maxwell’s Equations and Multisymplectic Scheme Supported by the Special Funds for Major State Basic Research Projects , G 1999, 032800

## 1 Introduction

Transient electromagnetic field problems are important problems in many modern technology applications. M. J. Gotay and J. E. Marsden $`etal`$ discussed the Lagrangian density for electromagnetism on a fixed background space-time $`X`$ with metric $`g`$, and discussed the existence of it’s Euler-Lagrange equation. In this paper we will discuss this problem from the point of inverse variational problem. As a complete depiction of the behavior of the electromagnetic field, Maxwell’s equations have very interesting structure. we will discuss its Lagrangian density and give its Multisymplectic Bridge’s(Hamilton) form.

The multisymplectic method is used to phrase the dynamics, which is described by a Euler-Lagrange system, in the terms of the $`finite`$-dimensional space of fields at a given event in space-time. Some Euler-Lagarange system can be formalized in a multisymplectic Hamilton form which manifest some characteristics of original systems, so we discretize multisymplectic Hamilton form directly in order to simulate original system intrinsically.

Determining whether a given nonlinear partial differential equation(PDE) is an Euler-Lagrange equation of some variational problem is fundamental work in finite-element method and nonlinear wave theory and etc.. In this paper we will use the straightforward result of this problem given by Atherton and Homsy to discuss self-adjointness of Maxwell’s equations in different representations. We show that the same equations in different forms shows different potentialness, consequently it can be expressed as different variational problems. So we choose an appropriate form of the Maxwell’s equations to derive it’s Lagrangian and then rewrite it in a multisymplectic Hamilton form by using the Legendre transformation.

This paper falls into 4 parts. In Sec. 2, we obtain the necessary and sufficient conditions of self-adjointness for a 1-st order form of the Maxwell’s equations in an inhomogeneous, isotropic and lossless medium. In Sec. 3, we rewrite the equations in a 2-ed order form by introducing two vector potential functions. A complex multisymplectic Hamilton form and a multisymplectic conservation law are given as well as a ordinary conservation law. In Sec. 4, as application of the multisymplectic method, a multisymplectic algorithm is constructed and numerical simulation is shown.

## 2 Lagrangian Formalisms for Maxwell’s Equations

First we shall introduce some basic concepts in variational problems. Let $`X^m`$, an open subset with smooth boundary $`X`$. $`M`$ is some function space. Given a vector differential equation with a vector of dependent variables $`u_i=u_i(x_1,x_2,\mathrm{},x_m)M,i=1,2\mathrm{},n`$, and a vector of differential operators $`N^j(u),j=1,2,\mathrm{},n`$, such that

$$N(u)=0.$$
(1)
We assume the Frech$`\stackrel{´}{e}`$t derivative of $`N`$ exist,

$$N_u^{}\varphi =\underset{h0}{lim}\frac{N(u+h\varphi )N(u)}{h}.$$
(2)
The arbitrary function $`\varphi `$ now belongs to $`M`$. The Frech$`\stackrel{´}{e}`$t derivative derives a matrix $`N_u^{}`$ with operator elements $`N^{}(i,j)`$ which are the derivatives of the operators $`N^i`$ with respect to $`u_j`$, so

$$N_u^{}=N^{}(i,j)_{i,j=1}^n.$$
(3)
It is called as differential operator matrix of operator matrix $`N`$. Then the self-adjoint test functional is given as

$$S(\psi ,\varphi )=_X\psi N_u^{}\varphi \text{d}v=_X\psi _iN^{}(i,j)\varphi _j\text{d}v.$$
(4)
From Vainberg’s theorem the condition of self-adjointness can be given as

$$S(\psi ,\varphi )=S(\varphi ,\psi )+_X\text{div}P\text{d}v,$$
(5)
if we set the boundary condition to be free, i.e. $`P`$ vanishes on $`X`$, then

$$N^{}(i,j)=N^{}(j,i)^{},$$
(6)
where $`N^{}(j,i)^{}`$ is denoted as the adjoint operator of $`N^{}(j,i)`$.

Surely a certain variational problem can lead to more than one Euler-Lagrange equation and different boundary conditions, responsibly the Lagrangian and the boundary condition are not unique in the inverse problems of finding corresponding functionals for a Euler-Lagrange equation, so the requirement for free boundary conditions in Ref. is too strict for some applications such as Multi-symplectic method which needs to verify self-adjointness of a given system to find one of its Lagrangians at first step of the method. For example, the wave equation such as

$$u_{tt}u_{xx}=0,$$
(7)
according to Ref. its self-adjointness condition is

$$\frac{N(u)}{u_2}|_X=0,\text{where}(t,x)X,u_2=u_{xx},u_{tt},$$
(8)
so it does not satisfy this condition whatever boundary condition is given, but it actually has Lagrangian, therefor we just set self-adjointness condition only to be condition (6) without free boundary condition.

Let us consider a $`k`$th order functional on maps $`u:^{1,3}^6`$, for Maxwell’s equation, $`u=[H_1,H_2,H_3,E_1,E_2,E_3]`$. Here $`^{1,3}`$ is a (1+3)-dimensional time-space domain, denoted by coordinates $`x_j`$, $`j=1,\mathrm{},4`$. $`^6`$ is the target space, denoted by coordinates $`y^i`$, $`i=1,\mathrm{},6`$, additionally we denote the space $`T_y^6T_x^{}^{1,3}`$ by $``$. The functional is defined by means of a Lagrangian $`L:`$ it is a function on the space $`^{1,3}\times ^6\times `$, i.e. $`L`$ is a smooth function defined on the bundle over $`^{1,3}\times ^6`$ with fiber over $`(x,y)`$ equal to $``$. We denote this bundle by $`T^6_{^{1,3}\times ^6}T^{}^{1,3}`$. Hence the functional is defined as

$$=_{^{1,3}}L(x_j,u^i(x),u_\nu ^i(x))\omega ,$$
(9)
where $`\omega `$ is some volume 4-form on $`^{1,3}`$, $`u_\nu ^i=\frac{u^i}{x_{j_1}\mathrm{}x_{j_\nu }}`$, $`\nu =1,2,\mathrm{},k`$. We define Legendre transformations as

$$\pi _i^\nu :=\frac{L}{u_\nu ^i}(x_j,u^i(x),u_\nu ^i(x)),$$
(10)
it leads to covariant Hamiltonian function $``$ on the multisymplectic manifold $`T^{}^6_{^{1,3}\times ^6}T^{1,3}`$, as

$$(x_j,u^i(x),\pi _i^\nu )=\pi _i^\nu u_\nu ^iL(x_j,u^i(x),u_\nu ^i(x))).$$
(11)
Now we consider the Maxwell’s equations in an inhomogeneous, isotropic, and lossless medium

$$\begin{array}{c}\hfill \mu \frac{H}{t}+\times E=K\\ \hfill \times H+\epsilon \frac{E}{t}=J,\end{array}$$
(12)
$`E=[E_1,E_2,E_3]^T`$ is the electric field strength, $`H=[H_1,H_2,H_3]^T`$ is the magnetic field strength with superscript T denoting transpose. $`\epsilon `$ is the permittivity, $`\mu `$ is the permeability, both $`\epsilon `$ and $`\mu `$ are scalar functions of space and time variables. $`J`$ is the external electric-current density, and $`K`$ is the external magnetic-current density. System (12) can be arranged in the form

$$𝒢𝒵=,$$
(13)
where $`𝒵=𝒵(x,y,z,t)`$ is the field vector consisting of the components of the electric field strength $`E`$ and the magnetic field strength $`H`$ as

$$𝒵=[H_1,H_2,H_3,E_1,E_2,E_3]^T$$
(14)
simplified into $`𝒵=[H,E]^T`$, similarly $`=(x,y,z,t)`$ is the source vector composed of the external electric current source $`J`$ and the external magnetic current source $`K`$. $`𝒢`$ is the operator matrix given by

$$𝒢=\left(\begin{array}{cc}\times & \epsilon _t\\ \mu _t& \times \end{array}\right),$$
(15)
where each element of $`𝒢`$ is a $`3rd`$ order matrix such as

$$\times =\left(\begin{array}{ccc}0& \frac{}{z}& \frac{}{y}\\ \frac{}{z}& 0& \frac{}{x}\\ \frac{}{y}& \frac{}{x}& 0\end{array}\right)\epsilon (\mu )_t=\epsilon (\mu )\left(\begin{array}{ccc}_t& 0& 0\\ 0& _t& 0\\ 0& 0& _t\end{array}\right).$$
(16)
Because the operator $`𝒢`$ is a linear operator matrix, i.e. each element of $`𝒢`$ is a linear operator, its differential operator matrix denoted as $`𝒢^{}`$ equals $`𝒢`$ itself. So the adjoint operator matrix of $`𝒢^{}`$ is as

$$𝒢^{}=\left(\begin{array}{cc}(\times )^\text{T}& \epsilon _t+\epsilon _t\\ \mu _t\mu _t& (\times )^\text{T}\end{array}\right),$$
(17)
It is obvious that the self-adjoint condition (6) is that $`𝒢^\text{T}=𝒢^{}`$, so the equations have to satisfy $`\epsilon =\mu `$, and the two coefficients are independent on time variable $`t`$. The system in vacuum of course satisfies $`𝒢^\text{T}=𝒢^{}`$. The systems with such self-adjointness condition have one Lagrangian as

$$\begin{array}{cc}\hfill L=& \frac{1}{2}<H,\times H>+\frac{1}{2}<E,\times E>\mu <H,E_t><H,J>\hfill \\ & +<E,K>,\hfill \end{array}$$
(18)
where $`\mu =\mu (x,y,z)=\epsilon (x,y,z)`$. How to calculate Lagrangian density form Euler-Lagrange equation is to be introduced in the next section.

Consider another form of operator matrix $`𝒢`$

$$𝒢_1=\left(\begin{array}{cc}\frac{1}{\epsilon }\times & _t\\ _t& \frac{1}{\mu }\times \end{array}\right),$$
(19)
with the source function $``$ changed in

$$_1=[\frac{1}{\epsilon }J,\frac{1}{\mu }K]^\text{T},$$
(20)
the dependent variable $`𝒵`$ is unchanged. But the differential operator matrix of $`𝒢_1`$ is not a self-adjoint operator matrix, since that

$$𝒢_1^{}=\left(\begin{array}{cc}A+\frac{1}{\epsilon }(\times )^\text{T}& _t\\ _t& B+\frac{1}{\mu }(\times )^\text{T}\end{array}\right)𝒢_1,$$
(21)
where $`A`$ and $`B`$ are two 3-rd order matrix as

$$A=\left(\begin{array}{ccc}0& (\frac{1}{\epsilon })_z& (\frac{1}{\epsilon })_y\\ (\frac{1}{\epsilon })_z& 0& (\frac{1}{\epsilon })_x\\ (\frac{1}{\epsilon })_y& (\frac{1}{\epsilon })_x& 0\end{array}\right),B=\left(\begin{array}{ccc}0& (\frac{1}{\mu })_z& (\frac{1}{\mu })_y\\ (\frac{1}{\mu })_z& 0& (\frac{1}{\mu })_x\\ (\frac{1}{\mu })_y& (\frac{1}{\mu })_x& 0\end{array}\right),$$
(22)
obviously if $`𝒢_1^\text{T}=𝒢_1^{}`$, $`\epsilon `$ and $`\mu `$ are free of space variables, additionally such form excludes the condition of $`\epsilon `$, $`\mu =0`$.

So we get another Lagrangian of system (12)

$$\begin{array}{cc}\hfill L=& \frac{1}{2\epsilon }<H,\times H>+\frac{1}{2\mu }<E,\times E><H,E_t>\frac{1}{\epsilon }<H,J>\hfill \\ & +\frac{1}{\mu }<E,K>.\hfill \end{array}$$
(23)
This Lagrangian seems equivalent to the Lagrangian (18), but in fact the above cases actually dealt with two different systems because of different restrictions for the coefficients $`\epsilon `$ and $`\mu `$ in their self-adjointness conditions.

Both the two Lagrangian are of first order functions, and the Legendre transformations are unnecessary in constructing their multisymplectic Hamilton form, because the conjugate momentums are $`H`$ and $`E`$ multiplying some coefficients, for example in $`L_2`$

$$\begin{array}{c}\hfill \frac{L_2}{\times H}=\frac{1}{2\epsilon }H,\\ \hfill \frac{L_2}{\times E}=\frac{1}{2\mu }E,\\ \hfill \frac{L_2}{E_t}=H.\end{array}$$
(24)
In fact, for linear 1-st order system which is self-adjoint, if its coefficients satisfy some conditions, then it can be written in multisymplectic Hamilton form directly.

## 3 Multisymplectic Hamilton Forms for Maxwell’s Equations

Introducing two vector functions $`U`$ and $`V`$ satisfying $`U_t=E`$ and $`V_t=H`$, so the system can be rewritten in form $`(2)`$ with different operator matrix $`𝒢`$, different dependent variable $`𝒵`$ and source function $``$

$$\begin{array}{c}\hfill 𝒢=\left(\begin{array}{cc}\mu _t^2& \times _t\\ \times _t& \epsilon _t^2\end{array}\right),𝒵=[V_1,V_2,V_3,U_1,U_2,U_3]^T,\\ \hfill =[K,J]^\text{T}.\end{array}$$
(25)
The differential operator matrix $`𝒢^{}`$ equals $`𝒢`$ itself. If $`\epsilon `$ and $`\mu `$ are not dependent on spatial variables, then it is easy to demonstrate that $`𝒢^{}`$ is self-adjoint, so the self-adjointness of the form (13) with sketch (25) is verified, and the potential functional is given by

$$\begin{array}{cc}\hfill F(V,U)& =_\mathrm{\Omega }_0^1[V,U]𝒢(\lambda 𝒵)\text{d}\lambda \text{d}v\hfill \\ & =_\mathrm{\Omega }(\frac{1}{2}\mu <V,V_{tt}>+\frac{1}{2}<V,\times U_t>+\frac{1}{2}\epsilon <U,U_{tt}>\hfill \\ & \frac{1}{2}<U,\times V_t>+UJ+VK)\text{d}v\hfill \\ & =_\mathrm{\Omega }(\frac{1}{2}\mu <V_t,V_t>\frac{1}{2}<V_t,\times U>\frac{1}{2}\epsilon <U_t,U_t>\hfill \\ & +\frac{1}{2}<U_t,\times V>+UJ+VK+\text{div}P)\text{d}v,\hfill \end{array}$$
(26)
where $`<,>`$ represents the inner production of vectors. We know the Lagrangian is unique for a given Euler-Lagrange equation under difference of a term $`\text{div}P`$. Here for simplicity we consider the Lagrangian as

$$\begin{array}{cc}\hfill L=& \frac{1}{2}\mu <V_t,V_t>+\frac{1}{2}<V_t,\times U>\hfill \\ & +\frac{1}{2}\epsilon <U_t,U_t>\frac{1}{2}<U_t,\times V>UJVK,\hfill \end{array}$$
(27)
correspondingly the generalized conjugate momentums are

$$\begin{array}{cc}\hfill P=\frac{L}{V_t}=\mu V_t+\frac{1}{2}\times U,& \frac{L}{\times V}=\frac{1}{2}U_t,\hfill \\ \hfill Q=\frac{L}{U_t}=\epsilon U_t\frac{1}{2}\times V,& \frac{L}{\times U}=\frac{1}{2}V_t.\hfill \end{array}$$
(28)
We get the covariant Hamiltonian by

$$\begin{array}{cc}\hfill S=& <P,V_t>+<Q,U_t>+<\frac{L}{\times V},\times V>\hfill \\ & +<\frac{L}{\times U},\times U>L\hfill \\ \hfill =& <P,H>+<Q,E>\frac{1}{2}\mu <H,H>\frac{1}{2}\epsilon <E,E>+UJ+VK,\hfill \end{array}$$
(29)
here $`Z=[H,E,V,U,P,Q]^T`$. So the Maxwell’s equations are transformed into the following form

$$\begin{array}{c}\hfill \frac{1}{2}\times U=P\mu H,\\ \hfill \frac{1}{2}\times V=Q\epsilon E,\\ \hfill P_t\frac{1}{2}\times E=K,\\ \hfill Q_t+\frac{1}{2}\times H=J,\\ \hfill V_t=H,\\ \hfill U_t=E.\end{array}$$
(30)
The above equations can be organized in the universal form

$$Z_t+𝒦\times Z=_ZS(Z),$$
(31)
where the rotation action $`\times Z`$ denotes $`[\times H,\times E,\times V,\times U,\times P,\times Q]^\text{T}`$, and $`𝒦\times Z`$ consists of three components such that

$$𝒦\times Z=𝒦_1Z_x+𝒦_2Z_y+𝒦_3Z_z.$$
(32)
So

$$=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& I& 0\\ 0& 0& 0& 0& 0& I\\ 0& 0& I& 0& 0& 0\\ 0& 0& 0& I& 0& 0\end{array}\right),𝒦_i=\left(\begin{array}{cccccc}0& 0& 0& \frac{1}{2}_i& 0& 0\\ 0& 0& \frac{1}{2}_i& 0& 0& 0\\ 0& \frac{1}{2}_i& 0& 0& 0& 0\\ \frac{1}{2}_i& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right),$$
(33)
where $`I`$ is the identity element belonging to $`R^{3\times 3}`$. $``$ is obviously anti-symmetric, and

$$R_1=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right)R_2=\left(\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right)R_3=\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right),$$
(34)
ultimately $`𝒦_i`$ belongs to $`R^{18\times 18}`$ and it is anti-symmetric. Additionally the rotation operator $`\times =R_1\frac{}{x}+R_2\frac{}{y}+R_3\frac{}{z}`$.

The representation (31) is simplified expression in vector form, its complete extension is the multisymplectic Hamiltonian system which was first introduced by Bridges and Derks in . It has the multisymplectic conservative law

$$\frac{}{t}𝝎+\times 𝜿=0,$$
(35)
$`𝝎`$ and $`𝜿`$ are the presymplectic forms

$$𝝎=\frac{1}{2}\text{d}Z\dot{}\text{d}Z,𝜿=\frac{1}{2}\text{d}Z\dot{}𝒦\text{d}Z,$$
(36)
where $`\text{d}Z=(\text{d}Z_1,\text{d}Z_2,\mathrm{})^\text{T}`$, and for two $`n`$-order vector functions $`p=(p_1,\mathrm{},p_n)^\text{T}`$ and $`q=(q_1,\mathrm{},q_n)^\text{T}`$, notation $`\dot{}`$ is defined as

$$\text{d}p\dot{}\text{d}q=\underset{i=1}{\overset{n}{}}\text{d}p_i\text{d}q_i.$$
(37)
The system (31) also has the local multisymplectic conservation law

$$\frac{1}{2}\times \left[\text{d}U\dot{}\text{d}H\text{d}V\dot{}\text{d}E\right]+\frac{}{t}(\text{d}V\dot{}\text{d}P+\text{d}U\dot{}\text{d}Q)=0.$$
(38)
Now we return to the second case introduced in Sec. 2. We can rewrite it in the multisymplectic Hamilton form directly

$$\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)\frac{}{t}\left(\begin{array}{ccc}& H& \\ & E& \end{array}\right)+\left(\begin{array}{cc}\frac{1}{\epsilon }_1& 0\\ 0& \frac{1}{\mu }_1\end{array}\right)\frac{}{x}\left(\begin{array}{ccc}& H& \\ & E& \end{array}\right)=\left(\begin{array}{ccc}& J/\epsilon & \\ & K/\mu & \end{array}\right),$$
(39)
here the covariant Hamiltonian is $`S=\frac{1}{\epsilon }<J,H>\frac{1}{\mu }<K,E>`$. It has the multisymplectic conservation law

$$\times \left[\frac{1}{\epsilon }\text{d}H\dot{}\text{d}H+\frac{1}{\mu }\text{d}E\dot{}\text{d}E\right]+\frac{}{t}(\text{d}E\dot{}\text{d}H)=0.$$
(40)
## 4 A Multisymplectic Scheme for the Maxwell’s Equations

Maxwell’s equations are simulated by many methods. Its first numerical simulation is given by Yee in 1966, which was developed in many more efficient algorithms. In this section, we adopt a multisymplectic scheme to approximate the system.

For simplicity, we consider $`1+1`$-dimensional multisymplectic form of (30), and suppose that the external magnetic-current vanished, so $`K=0`$, and the rotation operator becomes

$$\times Z=_1Z_x.$$
(41)
We discretize the equations (30) by using midpoint scheme in both $`t`$ and $`x`$ directions,

$$\begin{array}{cc}& _1\frac{U_{i+1j+\frac{1}{2}}U_{ij+\frac{1}{2}}}{2\mathrm{\Delta }x}=P_{i+\frac{1}{2}j+\frac{1}{2}}\mu _{i+\frac{1}{2}}H_{i+\frac{1}{2}j+\frac{1}{2}},\hfill \\ & _1\frac{V_{i+1j+\frac{1}{2}}V_{ij+\frac{1}{2}}}{2\mathrm{\Delta }x}=Q_{i+\frac{1}{2}j+\frac{1}{2}}\epsilon _{i+\frac{1}{2}}E_{i+\frac{1}{2}j+\frac{1}{2}},\hfill \\ & \frac{P_{i+\frac{1}{2}j+1}P_{i+\frac{1}{2}j}}{\mathrm{\Delta }t}_1\frac{E_{i+1j+\frac{1}{2}}E_{ij+\frac{1}{2}}}{2\mathrm{\Delta }x}=0,\hfill \\ & \frac{Q_{i+\frac{1}{2}j+1}Q_{i+\frac{1}{2}j}}{\mathrm{\Delta }t}+_1\frac{H_{i+1j+\frac{1}{2}}H_{ij+\frac{1}{2}}}{2\mathrm{\Delta }x}=J_{i+\frac{1}{2}j+\frac{1}{2}},\hfill \\ & \frac{V_{i+\frac{1}{2}j+1}V_{i+\frac{1}{2}j}}{\mathrm{\Delta }t}=H_{i+\frac{1}{2}j+\frac{1}{2}},\hfill \\ & \frac{U_{i+\frac{1}{2}j+1}U_{i+\frac{1}{2}j}}{\mathrm{\Delta }t}=E_{i+\frac{1}{2}j+\frac{1}{2}}.\hfill \end{array}$$
(42)
Here the subscript $`i`$ and $`j`$ denote respectively $`x`$ and $`t`$ directions. The above scheme is called the Preissman multisymplectic scheme which preserves the multisymplectic conservative law (35) in discrete space, such as

$$\frac{𝜿_{i+1j+\frac{1}{2}}𝜿_{ij+\frac{1}{2}}}{2\mathrm{\Delta }x}+\frac{𝝎_{i+\frac{1}{2}j+1}𝝎_{i+\frac{1}{2}j}}{\mathrm{\Delta }t}=0,$$
(43)
where

$$\begin{array}{cc}\hfill 𝜿_{ij+\frac{1}{2}}& =_1\text{d}U_{ij+\frac{1}{2}}\dot{}\text{d}H_{ij+\frac{1}{2}}_1\text{d}V_{ij+\frac{1}{2}}\dot{}\text{d}E_{ij+\frac{1}{2}},\hfill \\ \hfill 𝝎_{i+\frac{1}{2}j}& =\text{d}V_{i+\frac{1}{2}j+1}\dot{}\text{d}P_{i+\frac{1}{2}j+1}+\text{d}U_{i+\frac{1}{2}j+1}\dot{}\text{d}Q_{i+\frac{1}{2}j+1}.\hfill \end{array}$$
(44)
Eliminating the auxiliary variables $`V,U,P`$ and $`Q`$ from the Preissman scheme we get a nine-point multisymplectic integrator as

$$\begin{array}{cc}& _1\frac{E_{i+2j+2}+2E_{i+2j+1}+E_{i+2j}E_{ij+2}2E_{ij+1}E_{ij}}{4\mathrm{\Delta }x}\hfill \\ \hfill =& \mu _{i+3/2}\frac{H_{i+2j+2}+H_{i+1j+2}H_{i+2j}H_{i+1j}}{4\mathrm{\Delta }t}\hfill \\ & \mu _{i+1/2}\frac{H_{i+1j+2}+H_{ij+2}H_{i+1j}H_{ij}}{4\mathrm{\Delta }t},\hfill \\ & _1\frac{H_{i+2j+2}+2H_{i+2j+1}+H_{i+2j}H_{ij+2}2H_{ij+1}H_{ij}}{4\mathrm{\Delta }x}\hfill \\ \hfill =& \epsilon _{i+3/2}\frac{E_{i+2j+2}+E_{i+1j+2}E_{i+2j}E_{i+1j}}{4\mathrm{\Delta }t}\hfill \\ & +\epsilon _{i+1/2}\frac{E_{i+1j+2}+E_{ij+2}E_{i+1j}E_{ij}}{4\mathrm{\Delta }t}\hfill \\ & +(J_{i+3/2j+3/2}+J_{i+3/2j+1/2}+J_{i+1/2j+3/2}+J_{i+1/2j+1/2})/2.\hfill \end{array}$$
(45)
The numerical process can be depicted by the following expression

$$AZ^{j+2}=BZ^{j+1}+CZ^j+\overline{J},$$
(46)
where

$$A=\left(\begin{array}{ccccccc}D& A_0& 0& 0\mathrm{}& 0& 0& 0\\ B_0& D& A_1& 0\mathrm{}& 0& 0& 0\\ 0& B_1& D& A_2\mathrm{}& 0& 0& 0\\ \\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ \\ 0& 0& 0& 0\mathrm{}& B_{s1}& D& A_s\\ 0& 0& 0& 0\mathrm{}& 0& B_s& D\\ & \mathrm{}\end{array}\right),$$
(47)
$$C=\left(\begin{array}{ccccccc}D& A_0& 0& 0\mathrm{}& 0& 0& 0\\ B_0& D& A_1& 0\mathrm{}& 0& 0& 0\\ 0& B_1& D& A_2\mathrm{}& 0& 0& 0\\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ \\ 0& 0& 0& 0\mathrm{}& B_{s1}& D& A_s\\ 0& 0& 0& 0\mathrm{}& 0& B_s& D\\ & \mathrm{}\end{array}\right),$$
(48)
$$B=\left(\begin{array}{cccc}2D& & & \\ & 2D& & \\ \\ & & \mathrm{}& \\ \\ & & & 2D\end{array}\right),\overline{J}=\left(\begin{array}{ccc}& a(0)+b(0)& \\ & a(1)+b(1)& \\ & \mathrm{}& \end{array}\right),$$
(49)
where $`a(i)=J_{i+\frac{1}{2}j+\frac{1}{2}}+J_{i+\frac{3}{2}j+\frac{1}{2}},b(i)=J_{i+\frac{1}{2}j+\frac{3}{2}}+J_{i+\frac{3}{2}j+\frac{3}{2}}`$, and

$$D=\left(\begin{array}{ccccc}0& \mathrm{\Delta }t_1& & & \\ \mathrm{\Delta }t_1& 0& \mathrm{\Delta }t_1& & \\ & \mathrm{\Delta }t_1& 0& \mathrm{\Delta }t_1& \\ \\ & \mathrm{}& \mathrm{}& \mathrm{}& \\ \\ & & & \mathrm{\Delta }t_1& 0\end{array}\right).$$
(50)
Set $`r=1,2,\mathrm{}`$, then for every $`i=1,2,\mathrm{},s,\mathrm{}`$, there are

$$A_i=\left(\begin{array}{ccccccc}e_0+e_1& e_1& 0& 0\mathrm{}& 0& 0& 0\\ e_1& e_1+e_2& e_2& 0\mathrm{}& 0& 0& 0\\ 0& e_2& e_2+e_3& e_3\mathrm{}& 0& 0& 0\\ & & & \mathrm{}& \mathrm{}& \mathrm{}& \\ \\ 0& 0& 0& 0\mathrm{}& e_{r2}& e_{r2}+e_{r1}& e_{r1}\\ 0& 0& 0& 0\mathrm{}& 0& e_{r1}& e_{r1}+e_r\end{array}\right),$$
(51)
here $`e_j=\mathrm{\Delta }x\epsilon _{j+ir+\frac{1}{2}}IR^{3\times 3},j=1,2,\mathrm{},r`$. $`B_i`$ is defined similarly and just change $`e_j`$ in $`u_j=\mathrm{\Delta }x\mu _{j+ir+\frac{1}{2}}IR^{3\times 3}`$.

To check the integrator (45), we consider a simple example. We shall take $`\epsilon =\mu =1`$ and $`J=0`$. Furthermore, we take exact initial conditions as

$$\begin{array}{c}\hfill H_1(x,0)=H_3(x,0)=0,H_2(x,0)=\sqrt{\frac{\epsilon }{\mu }}\mathrm{sin}(x),\\ \hfill E_1(x,0)=E_2(x,0)=0,E_3(x,0)=\mathrm{sin}(x),\end{array}$$
(52)
under this initial condition we can get a exact solution as following

$$\begin{array}{c}\hfill H_1=H_3=0,H_2=\sqrt{\frac{\epsilon }{\mu }}\mathrm{sin}(x\sqrt{\frac{1}{\epsilon \mu }}t),\\ \hfill E_1=E_2=0,E_3=\mathrm{sin}(x\sqrt{\frac{1}{\epsilon \mu }}t),\end{array}$$
(53)
We compute $`H_2`$ in rectangles $`[0,2\pi +3]\times [0,0.1]`$, $`[0,2\pi +3]\times [0,1]`$ and choose $`\mathrm{\Delta }x=(2\pi +3)/61`$,$`\mathrm{\Delta }t=0.01`$. Fig.1 and Fig.2 shows the numerical solutions respectively. We can see the solution is moving to $`x`$-direction positively. In the rectangle $`[0,2\pi +3]\times [0,10]`$, the trend of the solution is more manifest.

The absolute error at t-number=1000 is shown in Fig. 4 by the ’error-line’, and the ’exact solution-line’ is the graph of solution $`0.01\times H_2`$. The former has a very regular route which is so like the solution graph, and their vertexes are almost at the zero points of each other. The regular similitude give a challenge to modify the scheme which is expected in further task.

For the multisymplectic form (39) there are many existing schemes to simulate the system and to preserve the multisymplectic conservation law (40), such as the usual midpoint scheme which doesn’t present here. For cases of more general $`\epsilon `$ and $`\mu `$, the discussion for their multisymplecticity is similar, so omitted here.


text:
# On the dynamical anomalies in the Hamiltonian Mean Field model

## I Introduction

The present effort is devoted to study the N-dependence in the thermodynamical variables and dynamical behavior of the well-known the Hamiltonian mean field model ant ; lat1 ; lat2 ; lat3 ; lat4 ; zanette . This model describes a system of *N* planar classical spins interacting through an infinite-range potential ant :

$$H_N(\theta ,L;I,g)=K_N(L;I)+V_N(\theta ;g)$$

$$\underset{i=1}{\overset{𝑁}{=}}\frac{1}{2I}L_i^2+\frac{1}{2}g\underset{i,j=1}{\overset{𝑁}{}}\left[1\mathrm{cos}\left(\theta _i\theta _j\right)\right]\text{,}$$
(1)
where $`\theta _i`$ is the *i-th* angle and $`L_i`$ the conjugate variable representing the angular momentum, $`I`$ is the moment of inertia of the rotator and $`g`$ the coupling constant. This is an inertial version of the ferromagnetic X-Y model stan , which interaction is not restricted to first neighbors but is extended to all couples of spins.

Traditionally the coupling constant $`g`$ is suppose to be *N*-dependent:

$$g=\frac{\gamma }{N},$$
(2)
assumption which makes $`H`$ formally extensive ($`VN`$ when $`N\mathrm{}`$), since the energy remains non-additive and the system can not be trivially divided in independent subsystems. A fundamental observable of this model system is the magnetization density $`𝐦`$, which is given by:

$$𝐦=\frac{1}{N}\underset{i=1}{\overset{𝑁}{}}𝐦_i\text{,}$$
(3)
where $`𝐦_k=[\mathrm{cos}\theta _k,\mathrm{sin}\theta _k]`$ is the intrinsic magnetization of the *k*-th particle.

The canonical analysis of the model predicts a second-order phase transition from a low-energy ferromagnetic phase with magnetization $`m1`$ to a high-energy paramagnetic phase with $`m0`$. The dependence of the energy density $`U=E/N`$ on the temperature $`T`$, the caloric curve, is given by:

$$U=\frac{T}{2}+\frac{1}{2}\left[1m^2\left(T\right)\right]$$
(4)
($`U`$ and $`T`$ in units of $`\gamma `$) and this is shown at the Fig.(1). The critical point is at the energy $`U_c=0.75`$ corresponds to the critical temperature of $`T_c=0.5`$ ant .

The dynamical behavior of the HMF has been extensively investigated. Microcanonical simulations of dynamics are in general in good agreement with the canonical solution, except for a region below $`U_c`$, where some anomalous phenomena have been observe, such as ensemble inequivalenceand the slow relaxation regimen towards the Boltzmann-Gibbs’ equilibrium lat2 ; lat3 ; lat4 ; zanette .

This ensemble inequivalence, characterized by the existence of a negative specific heat, appears as a dynamical feature. The origin of this anomalie is explained by the existence of certain *quasi-stationary state* lat3 ; lat4  which is not the canonical one and depends on the initial conditions in the microcanonical numerical computation of the dynamics. In this case, the system needs a long time to relax to the canonical equilibrium state because of the duration of this metastable state has a linear dependence with $`N`$. This fact has been interpreted by Tsallis as the non-commutativity of the thermodynamic limit ($`N\mathrm{})`$ with the infinite time limit ( $`t\mathrm{}`$): when the first is performed before the second, the system will not relax to the Boltzmann-Gibbs’ equilibrium.

Recent studies have shown that the slow-relaxation regime at long times is clearly revealed by numerical realizations of the model, but no traces of quasistationarity are found during the earlier stages of the evolution. This fact suggested the nonergodic properties of this system in the short-time range, which could make the standard statistical description unsuitable zanette .

However, it there could be a mistake in the understanding of the macroscopic description of this model. There are many studies in long-range Hamiltonian system where the validity of the thermodynamical description in the thermodynamic limit is intimately related to the scaling with $`N`$ of the thermodynamic variables and potentials. The problems with the nonextensive nature of those systems are avoided using the Kac prescription kac  in which the coupling constants are scaled by some power of $`N`$ in order to deal with an extensive total energy $`E`$. It is well-known that this kind of assumption affected the time scale of the system evolution, so that, its application should be taken with care.

## II Microcanonical analysis

M. Antoni, H. Hinrichsen and S. Ruffo carried out the microcanonical analysis of the HMF in the ref.ant2 . However, in their work they considered the equivalence of the microcanonical with the canonical ensemble in the thermodynamic limit, which depends on the validity of the canonical solution. In this section we will perform this analysis working directly on the microcanonical ensemble.

The microcanonical states density $`\mathrm{\Omega }`$ is calculated as follows:

$$\mathrm{\Omega }(E,N;I,g)=\frac{1}{N!}\frac{d^N\theta d^NL}{\left(2\pi \mathrm{}\right)^N}\delta \left[EH_N(\theta ,L;I,g)\right]\text{.}$$
(5)
The integration by $`d^NL`$ yields:

$$\frac{1}{N!\mathrm{\Gamma }\left(\frac{N}{2}\right)}\left(\frac{2\pi I}{\mathrm{}^2}\right)^{\frac{1}{2}N}\frac{d^N\theta }{\left(2\pi \right)^N}\left[EV_N(\theta ;g)\right]^{\frac{1}{2}N1}\text{.}$$
(6)
The microcanonical accessible volume $`W`$ is expressed as $`W=\mathrm{\Omega }g/2`$111For continuous variables, $`W`$ is only well-defined after a coarsed grained partition of phase space, which is the reason why it is considered an small energy constant $`\delta \epsilon `$ in order to make $`W`$ dimensionless. However, during the thermodynamic limit $`N\mathrm{}`$, the selection of $`\delta \epsilon `$ does not matter if it is small.. We are interested in describing the large $`N`$ limit, where $`W`$ is given by:

$$W\left(\frac{2\pi Ig}{N\mathrm{}^2}\right)^{\frac{1}{2}N}d^2𝐦\left[m^2+2U1\right]^{\frac{1}{2}N}f(𝐦;N),$$
(7)
where it was introduced the dimensionless parameter $`U`$:

$$E=gN^2U\text{,}$$
(8)
and the function $`f(𝐦;N)`$:

$$f(𝐦;N)=\frac{d^N\theta }{\left(2\pi \right)^N}\delta \left[𝐦\frac{1}{N}\underset{k=1}{\overset{𝑁}{}}𝐦_k\right]\text{.}$$
(9)
We estimate the function $`f(𝐦;N)`$ for large $`N`$ by using the *steepest decent method*. Expressing the delta function by mean of the Fourier’s integral representation, the Eq.(9) is rewritten as follows:

$$f(𝐦;N)=N^2\frac{d^2𝐤}{\left(2\pi \right)^2}\mathrm{exp}\left[iN\overline{𝐤}𝐦\right]\left[I_0\left(i\sqrt{\overline{𝐤}^2}\right)\right]^N$$

where $`\overline{𝐤}=𝐤+i𝐱`$, $`𝐱𝐑^2`$, and $`I_0\left(z\right)`$ is the modified Bessel function of zero order. It is easy to obtain the main contribution of the integral when $`N\mathrm{}`$ maximizing the integral function via the parameter $`𝐱`$. That is:

$$f(𝐦;N)N^2\mathrm{exp}\{N[xm\mathrm{ln}I_0\left(x\right)]$$

$$\frac{1}{2}\mathrm{ln}\left[\left(\frac{N}{2\pi }\right)^2\frac{m\left(x\right)}{x}\frac{d}{dx}m\left(x\right)\right]\},$$
(10)
where $`x=\left|𝐱\right|`$ is related with $`m=\left|𝐦\right|`$ by:

$$m=m\left(x\right)=\frac{I_1\left(x\right)}{I_0\left(x\right)}\text{.}$$
(11)
It is necessary to say that the validity of the steepest decent procedure is ensured by the positivity of the argument of the logarithmic term in the Eq.(10).

The main contribution to the entropy per particle $`s(U,N;I,g)=\mathrm{ln}W/N`$ for $`N`$ large can be obtained by using again the steepest decent method as follows:

$$s_o(U,N;I,g)=\frac{1}{2}\mathrm{ln}\left(\frac{2\pi Ig}{N}\right)+$$

$$+\underset{𝑥}{\mathrm{max}}\left\{\frac{1}{2}\mathrm{ln}\left[m^2\left(x\right)+\kappa \right]xm\left(x\right)+\mathrm{ln}I_0\left(x\right)\right\},$$
(12)
where $`\kappa =2U1`$. The maximization yields:

$$\left[\frac{m\left(x\right)}{x}m^2\left(x\right)\kappa \right]\frac{x}{m^2\left(x\right)+\kappa }\frac{d}{dx}m\left(x\right)=0\text{,}$$
(13)
which represents the energy dependence of magnetization density for the bulk matter:

$$\begin{array}{c}\begin{array}{ccc}\begin{array}{c}m=m(x)\\ U=\frac{1}{2}m\left(x\right)/x+\frac{1}{2}\left[1m^2\left(x\right)\right]\end{array}\}\hfill & \text{with}\hfill & x[0,\mathrm{})\hfill \end{array}\\ \text{and }m=0\text{ if }UU_c=0.75\end{array}.$$
(14)
The caloric curve (4) is obtained from the canonical parameter $`\beta `$:

$`\beta `$ $`={\displaystyle \frac{1}{T}}={\displaystyle \frac{s_0(U,N;I,g)}{U}}`$
$`=\{\begin{array}{ccc}x/m\left(x\right)\hfill & \text{with}\hfill & 0U<0.75\hfill \\ 1/(2U1)\hfill & \text{if}\hfill & U0.75\hfill \end{array}.`$ (17)
This is the same result found by M. Antoni, H. Hinrichsen and S. Ruffo in the ref.ant2 . Therefore, the subsequent analysis will be equivalent to their work. The microcanonical analysis confirms the existence of a second-order phase transition at the critical energy $`U_c=0.75`$ with a critical temperature of $`T_c=0.5`$.

It is interesting to note that although the main term of the entropy per particle is N-independent, there is an additional contribution which is N-dependent and makes the Boltzmann entropy nonextensive for an arbitrary N-dependence of the coupling constant $`g`$. Although the selection of the N-behavior of the coupling constant does not change the thermodynamical picture in the microcanonical ensemble, it affects the extensive or nonextensive character of the energy in the thermodynamic limit, as well as the temporal scale in the dynamical behavior of this model system (see in the next section).

In order to deal with an extensive Boltzmann entropy is necessary to assume the following thermodynamic limit for the Hamiltonian mean field model:

$$N\mathrm{},\text{ keeping constant }\frac{E}{N^3}\text{ and }\frac{g}{N}.$$
(18)
This assumtion defers from the one used by other investigators, where the N-dependence of the coupling constant chosen in order to deal with an extensive energy.

## III The origin of the slow relaxation regimen

According to the Hamiltonian of the model, its motion equations are given by:

$`\dot{\theta }_k`$ $`={\displaystyle \frac{1}{I}}L_k\text{,}`$
$`\dot{L}_k`$ $`=gN\left(m_y\mathrm{cos}\theta _km_x\mathrm{sin}\theta _k\right)\text{,}`$ (19)
where $`k=1,2,\mathrm{}N`$ and it was used the expression given in Eq.(3). It is convenient to introduce the following characteristic units for time and momentum:

$$\tau _0=\sqrt{\frac{I}{gN}}\text{ and }L_0=\sqrt{IgN}.$$
(20)
Note that, according to the thermodynamic limit (18), the characteristic time $`\tau _0`$ decreases with the $`N`$ increasing as $`\tau _01/N`$, while the characteristic momentum $`L_0`$ increases as $`L_0N`$. However, when it is assumed the *N*-dependence (2) for the coupling constant, the characteristic time $`\tau _0`$ is *N*-independent. The time scale of $`\tau _0`$ gives the characteristic time scale for the evolution of each rotator, so that, this time units is relevant to microscopic level. In the numerical computation of the dynamics this characteristic time scale should be taken into account in order to ensure a great accuracy of the energy conservation.

Using the above characteristic units, as well as introducing the tridimensional vectors $`𝐦_k=(\mathrm{cos}\theta _k,\mathrm{sin}\theta _k,0)`$ and $`𝐨_k=(0,0,L_k)`$, the motion equations (19) are rephrased as follows:

$$\dot{𝐦}_k=𝐨_k\times 𝐦_k\text{}\dot{𝐨}_k=𝐦_k\times 𝐦\text{.}$$
(21)
The magnetization density obeys the following dynamical equation:

$$\dot{𝐦}=𝚼=𝛀\times 𝐦+\mathrm{\Lambda }𝐦,$$
(22)
where $`𝚼`$ is the velocity of change of the magnetization density, $`𝛀`$, the angular velocity which characterizes the dynamical behavior of the $`𝐦`$ orientation, while $`\mathrm{\Lambda }`$ characterizes the behavior of its modulus. These quantities are given by:

$$𝚼=\frac{1}{N}\underset{k=1}{\overset{𝑁}{}}𝐨_k\times 𝐦_k,$$
(23)
$$𝛀=\frac{1}{Nm}\underset{k=1}{\overset{𝑁}{}}\left(𝐦_k𝐦\right)𝐨_k,\mathrm{\Lambda }=\frac{1}{Nm}\underset{k=1}{\overset{𝑁}{}}𝐨_k\times 𝐦_k𝐦.$$
(24)
Since $`𝐨_k\times 𝐦_k`$’s do not have a definite sign, the magnetization velocity of change $`𝚼`$ is a fast fluctuating quantity with a vanishing short time average (a time average along temporal interval comparable with the characteristic time of rotators evolution $`\tau _0`$) and standard deviation $`\sigma _𝚼`$ $`1/\sqrt{N}`$. This fact provokes that the short time average of the magnetization $`𝐦`$ varies slowly in comparison with the characteristic time scale of rotators evolution. Thus, there is two time scales in the dynamical evolution of the system: the microscopic time scale, $`\tau _{mic}=\tau _0`$, which is given by the characteristic time of the rotators evolution, and the macroscopic time scale, $`\tau _{mac}`$, which is given by the characteristic time of the short time average of magnetization density evolution, being this last considerably greater than the first time scale, $`\tau _{mac}\tau _{mic}`$.

We are interested in studying the dynamical evolution of the macroscopic observables for temporal scale comparable with $`\tau _{mac}`$. Since $`\tau _{mac}\tau _{mic}`$, at first approximation the dynamical evolution of $`𝐦`$ could be modeled by a Langevin equation with multiplicative noise celia ; foc :

$$\dot{𝐦}=𝚼=\sigma _1(m;U)\xi _1\left(t\right)𝐞_1+\sigma _2(m;U)\xi _2\left(t\right)𝐞_2.$$
(25)
whose parameters $`\sigma _1(m;U)`$ and $`\sigma _2(m;U)`$ could be estimated by using a microcanonical average with fixed magnetization density. This last estimation is a reasonable approximation since the magnetization density only exhibits very small fluctuations around its short time average for $`N`$large along time intervals comparable with $`\tau _{mic}`$. Here, $`𝐞_1`$ and $`𝐞_2`$ are unitary vectors in the parallel and perpendicular directions of $`𝐦`$ respectively, being $`𝐞_3=𝐞_1\times 𝐞_2(0,0,1)`$, the unitary vector along the third axis. $`\xi _1\left(t\right)`$ and $`\xi _2\left(t\right)`$ are Gaussian processes:

$$\xi _i\left(t\right)\xi _j\left(t^{}\right)=\tau _i\delta _{ij}\delta \left(tt^{}\right)$$
(26)
being $`\tau _i`$ the correlation times. The presence of $`\delta _{ij}`$ is due to the vanishing of the microcanonical average of $`𝚼_1𝚼_2=𝚼𝐞_1𝚼𝐞_2`$ with fixed magnetization density, which is easily proved by mean of a procedure analogue to the one used in the previous section for the determination of Boltzmann entropy. As already mentioned, the standard deviations of the components of the vector $`𝚼`$ depend on $`N`$ as $`\sigma _j(m;U)1/\sqrt{N}`$.

It is not difficult to understand that the correlation times $`\tau _i`$ are comparable with the microscopic time scale $`\tau _{mic}`$: the correlation times possess the same temporal scale of the characteristic time where the differenciable character of $`𝚼`$ is perceived, which is precisely characteristic time scale of rotators evolution.

According to the general theory, to this equation is associated the Fokker-Planck equation of the form:

$$\frac{}{t}F=_{m^\alpha }\left\{\frac{1}{2}\tau _j\sigma _je_j^\alpha _{m^\beta }\left[\sigma _je_j^\beta F\right]\right\},$$
(27)
which describes the dynamics of the magnetization density distribution function $`F=F(m,t;U)`$. In this expression, the $`e_j^\alpha `$’s represent the components of $`𝐞_j`$. Taking into account the N-dependence of $`\sigma _j`$, the right hand of the equation (27) decreases as $`1/N`$ during the $`N`$ increasing. Scaling the time scale as $`tt=Nt^{}`$, the Fokker-Planck equation becomes N-independent. This means that the characteristic time scale for the dynamical evolution of $`F(m,t;U)`$, as well as all those macroscopic observables derived from $`𝐦`$, have a linear dependence on $`N`$ during the increasing of system size (in $`\tau _0`$ units).

This is precisely the result which was observed in the numerical computation of the dynamics, where the relaxation time for the dynamical temperature $`T_D\left(t\right)=2K\left(t\right)/gN^2=2U1+m^2\left(t\right)`$ grows proportional to $`N`$ with the $`N`$ increasing. As already mentioned, the time scale of such dynamical behaviors is precisely the macroscopic time scale $`\tau _{mac}`$:

$$\tau _{mac}\tau _0N=\sqrt{\frac{IN}{g}}.$$
(28)
Thus, when the thermodynamil limit (18) is taken into account, the relaxation time is finite when $`N`$ is tended to infinity. Therefore, no anomalous dynamical behavior will be observed: the thermodynamic limit ($`N\mathrm{}`$) necesary for the ensemble equivalence between the microcanonical ensemble with the canonical one, will commute with the infinite time limit ($`t\mathrm{}`$) necesary for the equilibration of the temporal average of the physical observables. Thus, the slow relaxation regimen observed in the microcanonical numerical computation of dynamics appears as consequence of the consideration of the N-dependence of the coupling constant (2). The above analysis shows that the N-dependence of the coupling constant which ensures the extensivity of the Boltzmann entropy is more appropriate than the one which guarantees the extensivity of the energy. The first allows us to carry out a well-defined macroscopic description for the Hamiltonian mean field model, instead, the second leads to the existence of dynamical anomalies.

## IV Concluding remarks

As already showed, the microcanonical analysis revealed a thermodynamic limit (18) which defers from the traditional assumption of the $`N`$-dependence of the coupling constant (2). This thermodynamic limit is intimately related with the existence of the following self-similarity scaling behavior of the thermodynamical variables of the Hamiltonian mean field model:

$$\begin{array}{c}NN\left(\alpha \right)=\alpha N_,\\ EE\left(\alpha \right)=\alpha ^3E,\\ gg\left(\alpha \right)=\alpha g,\\ II\left(\alpha \right)=I,\end{array}\}W\left(\alpha \right)=[W\left(1\right),\alpha ]$$
(29)
where the functional $`[W,\alpha ]`$ defines an exponential self-similarity scaling laws for the microcanonical volume $`W`$:

$$[W,\alpha ]=\mathrm{exp}\left[\alpha \mathrm{ln}W\right],$$
(30)
being

$$W\left(\alpha \right)=W[E\left(\alpha \right),N\left(\alpha \right);g\left(\alpha \right),I\left(\alpha \right)],$$
(31)
where $`[W,\alpha ]`$ satisfies the self-similarity condition:

$$[[W,\alpha _1],\alpha _2]=[W,\alpha _1\alpha _2].$$
(32)
The reader may understand that this self-similarity scaling behavior is a generalization of the extensive properties exhibited by the traditional systems. As already discussed in our previous works, this self-similarity properties could be crucial in order to develop a well-defined thermodynamical description for a Hamiltonian nonextensive system vel1 ; vel2 ; vel3 . We showed in the ref. vel3  how self-similarity can be used in order to analyze the necessary conditions for the validity of Tsallis’ Nonextensive Statistics tsal1 . The results obtained in the present paper seem to be consistent with this idea too.

In future works we hope to obtain an analytical expression for the parameters appearing in the Fokker-Planck equation (27) in order to carry out a comparative study of this equation with the numerical computation of dynamics, since in this work we only concentrate on the analysis of their N-dependence.