text:
# Dark Matter in Galaxies: Conference Summary

## 1. Introduction

In the widely accepted Popperian interpretation of the scientific process, we proceed in two stages. First we use established theory and a mixture of observation, intuition and phantasy to set up a theory of how things work. Then we (or more likely our friends) make a determined effort to falsify our theory by finding a measurement or observation that is inconsistent with the theory. In most cases this effort succeeds fairly quickly, and we have to construct a new theory to continue the game. Occasionally resolute attempts to falsify the theory fail, and people gain confidence in the accuracy of its predictions. The theory is then considered established and becomes part of the infrastructure that is called upon in the first phase of the scientific process.

Sitting through this meeting, the conviction has grown on me that the Cold Dark Matter (CDM) theory has now reached the point at which it should be admitted as a Candidate Member, to the Academy of Established Theories, so that it can sit alongside the established theories of Maxwell, Einstein and Heisenberg and be used as a standard tool in the construction of new theories.

I start by summarizing the status of the CDM and MOND theories, followed by lists of headaches that the protagonists on each side have to confront, all interlarded with my current views on relevant issues. I finish with a ‘to do’ list. Since space is limited, I mostly rely on other contributors to the meeting for citations to individual papers: the word Bloggs is shorthand for “Bloggs, this volume and references therein.”

## 2. Cold Dark Matter’s Resumé

It is now about twenty years since cosmology became dominated by the theory that gravitational clustering is driven by CDM. This theory was honed by interaction with observations of galaxy clustering, the Ly$`\alpha `$ forest, the CMB and type Ia supernovae. At this meeting several speakers have concluded that astronomical data, especially those connected with gravitational lensing, agree with the predictions of the CDM theory. Recently the theory was rigorously tested by the WMAP satellite and passed brilliantly (Spergel). It will be further tested when data for additional years have been analysed. It is a theory with a significant hinterland in high-energy physics and as such can be tested in terrestrial experiments: attempts to detect the particles of cosmic CDM as they pass through the Earth are currently entering an exciting phase (Sellwood), and from around 2007 the LHC at CERN is expected to probe aspects of the underlying physics, which probably involves supersymmetry.

In the last few years it has been widely argued that there are conflicts between the predictions of CDM and observations of the internal structure of galaxies of various types. Much of this meeting has been taken up with discussion of these important questions. I argue below that these problems probably reflect difficulties in correctly deducing the predictions of the theory, rather than problems with the underlying physics.

## 3. Modified Newtonian Gravity’s Resumé

Milgrom made an important contribution to our subject by having the courage suggest that flat rotation curves reflected a failure of general relativity rather than the existence of vast quantities of unseen matter (Milgrom 1983). In a series of lucid and painstaking papers he explored in considerable depth this possibility, which he dubbed MOND. The essential ingredient of MOND is the assertion that there is a characteristic acceleration, $`a_0cH_0/2\pi `$, below which Newtonian gravity fails: when Newton predicts a gravitational acceleration $`g_\mathrm{N}<a_0`$, the actual acceleration is $`g_\mathrm{M}=\sqrt{g_\mathrm{N}a_0}`$.

Physicists are generally reluctant to give MOND the time of day because, despite the best efforts of Bekenstein & Milgrom (1984), a Lorentz covariant theory that has MOND for a Newtonian limit has not come to light. Astronomers tend to be more open-minded, and several studies have demonstrated that Milgrom’s MOND fitting formula enables one to predict the rotation curves of galaxies from their light profiles (Sancisi; Sellwood; Sanders & McGaugh 2002). What makes these fits impressive is that the mass-to-light ratios $`\mathrm{{\rm Y}}_C`$ that they assign to each galaxy vary with waveband $`C`$ and with galaxy type very much as one would expect on astrophysical grounds if the galaxies contained only stars and the observed gas (Bosma; de Jong). These studies have clearly established that a characteristic acceleration $`a_0`$ is involved in the relation between galaxy photometry and dynamics. This fact is at least as important for our understanding of galaxy structure and formation as the Tully-Fisher and Faber-Jackson relations.

There now seems to be little doubt that Einstein’s equations give an adequate account of the data only when a significant cosmological constant $`\mathrm{\Lambda }`$ appears in them. Physically, the non-zero value of $`\mathrm{\Lambda }`$ implies that the vacuum is a medium in which the energy density is non-zero, and $`\mathrm{\Lambda }`$ quantifies the density of this ‘dark energy’.111This is a poor name because dark energy comes in two forms. The simplest form is ‘dark matter’ (DM). This exerts negligible pressure. The form quantified by $`\mathrm{\Lambda }`$ is associated with a large negative pressure, i.e., a tension. So the physical quantity associated with $`\mathrm{\Lambda }`$ ought to be called ‘dark tension’ rather than dark energy. The existence of dark energy strongly suggests that our understanding of the way the vacuum responds to weak stimuli is seriously in error. Milgrom and others have suggested that there is a connection between the phenomenon of dark energy and the breakdown of the predictions of standard gravity when small accelerations are involved in the sense that both phenomena reflect differences between what actually happens and the predictions of Einstein’s theory with $`\mathrm{\Lambda }3(a_0/c)^2`$ set to zero (Milgrom 2002).

I find extremely seductive the idea that a better understanding of the nature of the vacuum would simultaneously explain the requirement for non-zero $`\mathrm{\Lambda }`$ and modify the way that particles moved in weak gravitational fields. MOND is a proposal for what the low-energy predictions of the correct dynamical model of the vacuum would look like. So while MOND itself could never become a member of the Academy of Established Theories, its parent theory most certainly could.

Thus CDM and MOND do not compete on equal terms: the former is a natural outgrowth of established physics, whilst the latter is a thing apart, a single visible peak of a mountain range that is otherwise enshrouded in cloud; a peak that is connected to the known world in an unkown way. The result is that in our efforts to understand how the Universe got to its present state, CDM can be used to form hypotheses and make predictions in a way that MOND cannot.

## 4. Headaches for MOND

### 4.1. Falling circular-speed curves

We heard that the Planetary Nebula Spectrograph on the William Herschel telescope has measured velocity dispersions in the outer parts $`R>2R_\mathrm{e}`$ of three intermediate-luminosity elliptical galaxies and find them all to fall in Keplerian fashion (Romanowsky). These data do not absolutely require the circular speeds of these galaxies to be Keplerian, but this is the most plausible interpretation of the data. Only two other intermediate-luminosity elliptical galaxies have had their stellar kinematics probed at $`R>2R_\mathrm{e}`$ and only one of these (NGC 2434; Rix et al., 1997) shows evidence for DM. Milgrom & Sanders (2003) point out that these very high surface-brightness galaxies generate accelerations around $`R_\mathrm{e}`$ that are an unusually large multiple of $`a_0`$. Consequently, MOND is expected to push the the circular-speed curve above Newton’s prediction only after a significant section of Keplerian fall. Hence, the work reported by Romanowsky is perhaps a triumph for MOND.

### 4.2. Non-universality of $`a_0`$

The value of $`a_0`$ required to banish DM from rich clusters of galaxies is about a factor of two larger than that fitted to data for disk galaxies (Sanders & McGaugh 2002 and references therein). Gerhard argued that the highest-quality data for elliptical galaxies at $`R<2R_\mathrm{e}`$ implies that DM manifests itself at accelerations $`a10^9\mathrm{m}\mathrm{s}^2`$ that are about an order of magnitude larger than the value of $`a_0`$ that is inferred from the dynamics of disk galaxies (Gerhard et al. 2001). However, this conclusion would appear to conflict with the subsequent results from planetary nebulae, which demonstrate that in two of the galaxies in the Gerhard et al. sample (NGC 3379 and NGC 4494), even 4 to 6 times further out than Gerhard et al. were able to go, there is no convincing evidence of DM. Baes & Dejonghe (2002) point out that scattering by dust at large radii of photons emitted near the centre can generate the kinematic signature of a DM halo when none is really present. Perhaps this effect is significant at the largest radii reached by Gerhard et al.

### 4.3. Halo flattening

Weak lensing studies, like studies of polar rings, suggest that the gravitational potentials of disk galaxies are significantly flattened. Is this a problem for MOND, as Hoekstra argued? It was not clear to me that those who answer ‘yes’ have borne in mind the strange behaviour of multipoles in MOND. Milgrom (1986) shows (i) that you cannot determine the quadrupole of the potential from the only quadrupole of the mass distribution, and (ii) that outside the body the quadrupole of the potential decays as $`r^\sqrt{3}`$, rather than $`r^3`$ as in Newton’s theory. Moreover, weak lensing is generated by accelerations that are only a few percent of $`a_0`$. To impart the velocity $`600\mathrm{km}\mathrm{s}^1`$ of the Local Group with respect to the CMB in a Hubble time, takes an acceleration $`acH_0/500=0.012a_0`$. Hence a single galaxy is unlikely to dominate the deflections probed by weak lensing, and non-circular isopotentials are not unexpected in MOND.

### 4.4. Merging galaxies

A famous photograph by Schweizer (1982) left little doubt that the merger of two disk galaxies of comparable mass yields an elliptical galaxy. The photograph shows the two long tidal tails of NGC 7252, together with the galaxy’s nearly relaxed core. The timescale on which the tidal tails evolve will be the same regardless of whether galactic rotation curves are kept flat by massive dark halos or modification of the law of gravity, because speed and distance are both fixed by the observations. Thus the time since the galaxies came into close contact is known. Schweizer showed the brightness distribution of the core obeys the $`R^{1/4}`$ law that is characteristic of elliptical galaxies. Thus the nuclei of the two galaxies have already completely merged. Simulations show that the nuclei can only spiral together in the time available if they can effectively surrender their energy and angular momentum to dark halos (Carignan; Barnes 1988). If we banish the halos by modifying the law of gravity, the galactic nuclei take much longer to merge because the vacuum cannot relieve them of their energy and angular momentum.

## 5. CDM and Galaxy Formation

History shows that even when we have the correct theory, prediction can be very difficult, even problematic: Newton published the Principia in 1686 but it was not until 1882 that Newcomb gave the now accepted value for the Newtonian precession of the perihelion of Mercury, and it is only ten years ago that Wisdom and others showed that there is significant chaos in the solar system; Euler wrote down his equation for hydrodynamics in 1755, but even in the sixth edition of his magisterial survey Hydrodynamics Lamb (1932) displayed a very incomplete understanding of shock fronts and in Art 284 ridiculed the Hugoniot jump condition. It is only in the last quarter century, with the development of chaos theory, that we have begun to understand the devastating impact that the non-linearity of Euler’s equations has on the equations’ predictive power.

CDM is a theory of the invisible. Apart from experiments designed to detect wimps in the laboratory and gravitational lensing work, tests of the CDM theory centre on the effect that CDM has on baryons. Hence galaxy-scale tests of CDM are inextricable entwined with the theory of galaxy formation, which nobody imagines is well understood. Because baryons interact with each other and electromagnetic radiation in tremendously complex ways, much cosmology is done with just CDM and dark energy. We do however now think that the mean density of baryons is $`\frac{1}{5}`$ that of DM, so the exclusion of baryons from the calculations is far from safe even on large scales. And baryons contribute much more heavily to the density precisely at the locations, in or near galaxies, where we make observations. So CDM-only simulations are, prima facie, unlikely to be reliable guides to the configuration of DM in and around galaxies. It’s not going to be easy to wring reliable predictions for galaxy-scale phenomena from the CDM model. Just as the fact that it is difficult to calculate the long-term dynamics of the solar system is no criticism of Newtonian dynamics, so the difficulty of predicting the DM distribution expected in the Milky Way is no criticism of the CDM theory.

Pure CDM simulations and galaxy-formation models tend to be done by the same workers in the same computers, but they are sharply different activities. The dynamics of pure CDM is reasonably clear-cut and well understood. We can feel confidence in the predictions it yields for large-scale phenomena, to which baryons contribute only modestly. These predictions are in excellent agreement with observation (Spergel). I have the impression that the predictions of pure CDM theory are in conflict with observation only where baryons either are, or likely have been, dynamically important. During this meeting I have come to the conclusion that we should treat these conflicts not as falsifications of CDM theory, but as pointers to failings in our understanding of galaxy formation. We are not even sure how baryons should be accumulated at the centres of DM halos in the purely dissipative case: the widely adopted prescription of conserving the adiabatic invariants of the DM as the potential gets deeper is probably quite wrong in the likely case that baryons accumulate through a series of mergers (Primack).

There is abundant evidence that non-gravitational energy is important for the dynamics of baryons that fall into DM halos of a wide range in mass, from dwarf galaxies right up to groups and even clusters of galaxies. The dynamics associated with this non-gravitational energy is currently no more than speculation.

## 6. Headaches for CDM

Many of the difficulties that CDM faces on small scales derive from the result that the centres of pure CDM halos are very cuspy. Even in the case of pure CDM, I believe we don’t understand properly how cores form, or what the smallest radius is at which we can trust the simulations. There are two problems: (i) the simulations still don’t reach as far down the mass function as the observations do, and (ii) we have an inadequate understanding of the role of discreteness effects. In the simulations the latter are important at early times when the Zel’dovich waves break, and they are distinct from the effects of two-body relaxation.222I used to worry about the latter, but I feel that Binney & Knebe (2002) showed that two-body relaxation is unimportant – see however Diemand et al. (2003) for a different point of view. However, through the cries and smoke of battle, something like a consensus seems to be emerging from the groups that do large-scale simulations of gravitational clustering (Navarro): in a pure-CDM halo the density profile rises in to the smallest resolvable radii with a slope $`\alpha =\mathrm{d}\mathrm{log}\rho /\mathrm{d}\mathrm{ln}r`$ that gradually lessens inwards, but is still $`>1`$ when last believable. The density profile of an individual halo depends weakly or not at all on the power spectrum from which the simulation starts (Knebe et al. 2002).

### 6.1. LSB galaxies

Stimulated by this prediction, observers have invested significant effort in determining the density profiles at the centres of the most DM dominated galaxies, since these clearly provide the cleanest tests. Although the controversy has yet to die, my impression is that they have shown that Low Surface Brightness (LSB) galaxies are unlikely to have such cuspy central mass profiles as the pure CDM simulations predict (Bolatto, Bosma, de Blok, Gentile, Mateo, Swaters, Trott). Dwarf spheroidal galaxies show the highest degree of DM domination (Wilkinson), and the conclusion of Kleyna et al. (2003) that the UMi dwarf has a nearly harmonic core is certainly a headache for CDM, but Draco does have a cuspy potential (Kleyna et al. 2002). It seems that in many, perhaps all LSB galaxies, the potentials are less cuspy than pure-CDM theory predicts. Will this problem be resolved at some future date by correctly adding baryons to the CDM simulations? The resolution will be hardest for the dSph galaxies, because their present baryon content is negligible. However, these are precisely the galaxies that are expected to lose most gas during galaxy formation, so it is not clear that their baryon contents were always negligible.

### 6.2. Halo substructure

In the last few years the worry has been abroad that DM halos are predicted to contain more substructure than the observed number of satellite galaxies suggests exists. It seems that this problem may have gone away: the number of objects predicted is a steep function of the sub-halo’s peak circular speed. Mapping this circular speed onto quantities observed for satellites, such as luminosity or central velocity dispersion, is non-trivial. It is now argued that when this mapping is done correctly, the CDM model correctly predicts the observed number of satellite galaxies (Primack). Moreover, in one interpretation of phenomena associated with strong gravitational lensing, a high level of halo substructure is inferred (Schneider, Mao and below).

### 6.3. HSB galaxies

Many arguments now indicate that the centres of high surface brightness (HSB) galaxies are baryon-dominated:

* Gerhard et al. (2001) find that the dynamics of elliptical galaxies is well explained by assuming that DM makes a negligible contribution to their central mass densities. The mass-to-light ratios then required are in the range expected on astrophysical grounds and vary with the colour of the galaxy in the expected way (Gerhard). Significant contributions to the mass density from DM would spoil this picture.
* The shapes of rotation curves are intimately connected to the underlying stellar light profile (Sancisi). The rotation curves of early-type disk galaxies peak at extremely small radii and must be dominated by stellar mass over most of the radial range probed (Noordermeer).
* In disk galaxies we more often than not see a fast bar (one that extends to $`0.8`$ of its corotation radius). The timescale for this to surrender most of its angular momentum to a slowly-rotating, embedding dark halo is $`f\frac{1}{2}(x+1/x)`$ times the dynamical time, where $`x=\rho _{\mathrm{bar}}/\rho _{\mathrm{DM}}`$ (Athanassoula). Since the high frequency of bars at the centres of galaxies implies that $`f`$ is quite large, $`x`$ must be far from unity.
* The non-axisymmetric velocities that a bar drives in surrounding gas are proportional to the mass of the bar, while the mass of the halo is limited by the overall rotation curve and the masses of bar and disk. Near maximal disks are required to generate velocities as high as those observed in some external galaxies (Bosma, Weiner). Similarly, the Milky Way’s ‘forbidden’ velocity features, such as the $`3`$kpc arm in the $`(l,v)`$ plots for CO and HI, require all available mass in the inner few kiloparsecs to be associated with the bar and the disk (Bissantz, Englmaier & Gerhard, 2003; Fux, 1999).
* We can weigh the Milky Way’s disk near the Sun and all recent investigators conclude that its entire surface density of $`41\mathrm{M}_{}\mathrm{pc}^2`$ can be accounted by known stars and interstellar gas (van Altena; Crézé et al. 1998; Holmberg & Flynn 2000; Olling & Merrifield 2001). Given that the total column density within $`1.1`$kpc of the plane is only $`71\mathrm{M}_{}\mathrm{pc}^2`$ (Kuijken & Gilmore, 1991), one can show that even as far out as the solar neighbourhood ($`R_0/R_\mathrm{d}>2.7`$), DM makes a smaller contribution than stars to the local circular speed (Binney & Evans, 2001).
* The microlensing optical depth to the Galactic centre has to be due to stars. Even the much reduced current values ($`\tau _610^6\tau 1.5`$: Popowski et al. 2001; Alfonso et al. 2003) can be explained only if all the mass that the circular-speed curve can accommodate in the central few kiloparsecs is invested in stars (Bissantz & Gerhard 2002).

I’d like to abuse my privilege of having the last word to contribute to the controversy about the implications of measurements of $`\tau _6`$. I coauthored two papers (Binney, Bissantz & Gerhard 2000; Binney & Evans 2001) that used an upper limit on the optical depth that can be achieved with a given mass when it is distributed smoothly around an ellipse. In the first paper we argued that the values, $`\tau _63`$, that were then current were physically impossible. Obligingly, Popowski et al. (2001) shortly afterwards reduced to $`2\pm 0.4`$ the most reliable value of $`\tau _6`$, that for red clump stars, which are bright enough for blending not to be a problem. In the second paper we showed that even the optical depth of Popowski et al. places a tight constraint on the power index $`\alpha `$ of the dark halo’s density profile between the Sun and the centre because we know the local density of DM, and $`\tau _6`$ limits the density of DM at the centre.

Recently the EROS collaboration have reported $`\tau _6=1.08\pm 0.3`$ (Alfonso et al. 2003) and at this meeting Merrifield has argued that flattening the dark halo to axis ratio $`q=0.8`$ allows $`\alpha `$ to be as large as unity without lowering $`\tau _6`$ below its likely value. Hence I hear it said on all sides that the Galactic microlensing data are now compatible with an NFW halo.

This conclusion is very wrong! What has been lost in the debate is how absurd the distributions of stars considered in our two papers are. These distributions are the ones that maximize the optical depth along $`b3.8^{}`$ for a given mass subject only to the constraints (i) that the density decreases exponentially with distance $`z`$ from the plane, with a scale-height $`z_0(R)`$, (ii) that the matter distribution is elliptical, with unconstrained ellipticity and a principal axis inclined at $`20^{}`$ to the Sun-centre line, and (iii) that the disk’s radial density profile is exponential. The resulting scale heights $`z_0`$ increase linearly with distance $`D`$ from the Sun so that $`z_0/D=|b|`$, the latitude at which the value of $`\tau _6`$ is set. The Galaxy’s stars are not distributed in this way! We used these absurd models to make the point that with the old data there was no way in which the lenses could be distributed through the Galaxy to achieve the required optical depth. Now that the microlensing measurements have been refined and are producing plausible values, it is time to add the constraint that the lenses, which are almost certainly ordinary stars, are distributed like the starlight. Since this distribution does not maximize the optical depth along a particular line of sight from the Sun, it produces a smaller optical depth per unit stellar mass than Binney & Evans assumed. The analysis of the COBE near-IR light distribution by Bissantz & Gerhard (2002) provides our best estimate of the distribution of stars near the centre. If these stars contain all the mass that is allowed in the centre, Bissantz & Gerhard conclude that $`\tau _6=1.27`$. So at the present time there is very little room for DM in the inner few kiloparsecs, even with the new lower optical depths.

## 7. Things to do

I have argued that ab initio prediction of the DM distribution in galaxies is hard because it involves solution of the full galaxy-formation problem. If we accept this conclusion, we should concentrate on mapping DM in galaxies, in an open frame of mind. This enterprise has two arms: (i) use of all available dynamical probes (stars, X-ray gas, atomic and molecular gas, gravitational micro- and macro-lensing) to map total mass, while (ii) using the most sophisticated models of stellar populations and interstellar gas to infer baryonic mass. As we have heard from many speakers, arm (i) is rather fully developed. By contrast, arm (ii) requires much more work, especially as regards the ISM. I feel the work described by de Jong is a valuable step along the right path.

We must go after cold molecular gas like bloodhounds, while avoiding a false dichotomy between ‘missing mass is all in DM’, and ‘missing mass is all in the ISM’. In the LSB regions of galaxies it is natural that gas should condense into very cold clouds that are hard to see in emission (Pfenniger). Observations in the $`21`$cm line have long told us that galaxies have HI disks that extend into extremely LSB regions. If the ISM in these regions were several times more massive than its HI component, we could understand a number of things:

* The existence of spiral structure in the extended disk around NGC 2915 (Masset).
* Why polar-rings lie above the Tully-Fisher relation for normal spirals, while their embedded disk galaxies lie on it (Arnaboldi).
* The existence of the baryonic Tully-Fisher relation (Combes; Matthews et al. 1998; McGaugh et al. 2000).
* The presence $`20`$kpc from the centre of M31 of a surprising amount of dust and star formation. If stars are forming, there must be unobserved H2 present, and if the ISM has the low metallicity expected so far out, the observed level of reddening is hard to understand without a significant column of H2 (Allen, Pfenniger).

Data from the new X-ray observatories have yet to have a big impact on studies of individual galaxies because they require extremely long exposure times, and we need to be more sophisticated in our modelling of hot gas. In particular, models need to include the spin of the body of trapped gas. I think there is considerable scope for linking our knowledge of stellar populations to the composition and global dynamics of the ISM.

We need to clarify the situation regarding DM in elliptical galaxies. Studies of gravitational lensing suggest that at $`2R_\mathrm{e}`$ both stars and DM make significant contributions (Schneider, Schechter). The flux ratios in the images of strongly lensed background objects require the mass distribution in the lensing galaxy to be lumpy (Mao). Schneider took the lumpiness to be sub-structure in the DM halo, such as that predicted by the DM clustering simulations. Schechter argued that the lumpiness of the stellar contribution to the mass distribution was responsible to the observed flux ratios. Observations of the variability of narrow and broad emission lines in the background source will resolve this issue within a few years and greatly clarify the distribution of DM in elliptical galaxies.

If we exclude cluster-centre galaxies, the dynamics of stars provides at best weak evidence for DM: the data favour some DM at $`R2R_\mathrm{e}`$, but what there is does not prevent the circular speed dropping in near Keplerian fashion at $`R>2R_\mathrm{e}`$ (Gerhard, Romanowsky). Compact DM halos that yield falling rotation curves have long been favoured by studies of tidal-tail formation and evolution (Dubinski et al. 1999 & references therein). There must be a worry that the rather different picture provided by the lensing data is biased: even if only a minority of ellipticals have massive dark halos, nearly all the observed lenses will belong to that minority.

Last but by no means least we have to press on with developing our understanding of how galaxies form. This is going to be a long job, but an immensely worthwhile one. We’ll probably crack it soonest if we accept CDM as background theory.


format: text

text:
# Using Artificial Intelligence for Model Selection

## 1 Introduction

In this paper we apply a new method for model selection on large data sets using the artificial intelligence algorithm Adaptive Simulated Annealing (ASA). We focus here on an epidemiological setting, although our techniques could be applied more widely. The objective is to analyze a large amount of data on many different characteristics of a given population and to select the model that best predicts the probability that an individual with various traits will have a particular disease outcome. The inclusion of categorical predictors that split up into multiple “dummy” variables and the possibility of cross terms between the different characteristics make the number of variables in a complete model prohibitively large. Thus it is necessary to select a small number of variables that will give the most informative model. Traditionally, this is done using forward regression or backward regression. Instead, we place an upper bound on the number of variables we want in our final model in advance and use ASA to decide which characteristics (or combinations of characteristics) each variable should represent to produce the model with the lowest $`C_p`$ statistic, a reflection of bias and variance. We ran computer simulations on populations of 100,000, 500,000, and 1,000,000, equivalent to the data that might be available for a major metropolitan area. Our results suggest that this method produces a model with a $`C_p`$ statistic that is consistently close to optimal whereas both forward and backward regression occasionally do not.

## 2 The Problem

We assume that we have access to a large data set with complete records of a population with a variety of fields, including both continuous variables (such as age and alcoholic consumption) and categorical variables (such as gender, ethnicity, and blood type). We also assume that we know whether each individual in the population has a particular disease, for example, liver cancer. (This is a binary response, since the outcome is either “yes” or “no”, but our method could be easily modified to predict a continuous response, for example, years of life lost due to liver cancer.) In our computer simulations we used models of varying numbers of people with characteristics randomly generated according to standard probabilities. We assigned probabilities of getting the disease according to various risk factors among the characteristics and let the computer randomly choose which individuals to afflict. The characteristics we used were age, gender, Hepatitis B viral infection, Hepatitis C viral infection, aflotoxin exposure, genetic marker, alcohol, and tobacco. Simulations with 1,000,000 people, equivalent to a major metropolitan area, resulted in about 400-450 people getting the disease in any given trial. We then selected an equal number of healthy people randomly from the population and added these to the diseased people to form a data set for a simulated population-based case-control study.

A classical problem of public health is then to develop a model from the data set that predicts an individual’s probability of getting the disease based on his/her combination of characteristics. To do this, we create a set of variables $`\stackrel{}{w}=\{w_i\}`$ as follows: For each continuous characteristic (such as age or alcohol consumption) we create one variable, and for each categorical characteristic (such as gender or blood type) with $`n`$ categories, we create $`n1`$ indicator variables that can take values of 0 or 1. The $`i`$th indicator formed from a particular categorical variable will indicate inclusion in category $`i`$. (We only need $`n1`$ of these indicator variables because if all such variables are 0, then we are guaranteed inclusion into the final category.) We then posit a model of the form

$$g(\stackrel{}{w})=\gamma _0+\underset{i}{}\gamma _iw_i+\underset{i<j}{}\gamma _{i,j}w_iw_j$$

in which the $`\gamma _i`$ and $`\gamma _{i,j}`$ are constant coefficients and we omit the cross term $`w_iw_j`$ when $`w_i`$ and $`w_j`$ arise from the same categorical variable. This model is linear (in the sense that if all variables but one are held constant, it is linear in the nonconstant variable), but we could as easily include higher degree terms and cross terms involving more than two variables if necessary. In practice it is uncommon for biological models to include terms involving three or more variables; we exclude these terms because they are difficult to interpret biologically.

To simplify the notation, we can change variables from $`\stackrel{}{w}`$ to $`\stackrel{}{x}`$, where each $`x_i`$ represents either 1 (the constant term), one of the $`w_i`$’s, or one of the $`w_iw_j`$’s. Then, replacing the various $`\gamma `$’s with a single set of coefficients $`\{\beta _i\}`$, the model above becomes

$$g(\stackrel{}{w})=f(\stackrel{}{x})=\underset{i}{}\beta _ix_i=\stackrel{}{\beta }\stackrel{}{x}$$
(1)
In the traditional approach to this problem, which, as we will see below, uses forward or backward regression, the ordering of the $`x_i`$’s, in other words, which $`x_i`$ is assigned to which $`w_i`$ or $`w_iw_j`$, is not important. However, our method using Adaptive Simulated Annealing depends heavily on finding a structured way of assigning the $`x_i`$’s so that $`x_i`$ will have some relation to $`x_j`$ when $`i`$ is close to $`j`$. We will discuss this further in Section 4 below.

One final modification is necessary to our model. To account for the fact that the response from a known individual is either “healthy” or “sick”, we can, at the end, transform the model (1) by the logistic function

$$P(\stackrel{}{x})=\frac{1}{1+e^{f(\stackrel{}{x})}}.$$

For a detailed explanation of this transformation, see Section 12.12 in .

The problem now is that of economy of terms; in our simulation, for example, the eight characteristics we studied produced 61 variables, $`x_0`$ to $`x_{60}`$. We therefore try to find a model using fewer terms and identify the $`x_i`$’s (and hence the individual characteristics or combinations thereof) that produce the most informative results. A commonly used measure of the efficiency of a model is the $`C_p`$ statistic

$$C_p=p+\frac{(s^2\widehat{\sigma }^2)(np)}{\widehat{\sigma }^2},$$

where $`p`$ is the number of variables in the current model, $`s^2`$ is the mean square error, and $`n`$ is the total possible number of variables. Finally, $`\widehat{\sigma }^2`$ is an estimate of $`\sigma ^2`$, the error variance in the population; since this latter quantity is unavailable, we take $`\widehat{\sigma }^2`$ to be the mean square error for the most complete model, i.e. the one incorporating all the variables. Then a lower $`C_p`$ value indicates a more desirable model. For a general overview of these terms, see Chapters 11 and 12 in ; for a more detailed development of the the $`C_p`$ statistic, see .

Two traditional methods of solving this problem are forward regression and backward regression. In the former, one starts with only the constant variable and adds variables to the model one at a time, each time selecting the variable that lowers the $`C_p`$ statistic the most. The process terminates when no variables can be added that lower the statistic. Backward regression is similar, starting with all the variables and progressively deleting them. Both procedures are computationally intensive, since they involve at each step computing a least squares solution to $`\text{X}\stackrel{}{\beta }=\stackrel{}{y}`$ (where X represents the data matrix) for each variable that is under consideration to be added or deleted from the model. This in turn requires finding the inverse (or pseudoinverse) of the matrix $`\text{X}^T\text{X}`$, which has dimension equal to the number of variables currently under consideration (up to 61x61 in our simulations). As greedy algorithms, they are also not guaranteed to produce an optimal solution.

As we noted above, since forward and backward regression test at each step every remaining variable to see which one should be added or deleted, these methods are not dependent on the ordering of the variables $`x_i`$.

Our new method described in Section 4 below is to determine an appropriate number of variables in advance and then to use ASA to determine which variables to include in the model to yield the lowest $`C_p`$ statistic. We will see that the ordering of the $`x_i`$’s is important to this method.

## 3 Adaptive Simulated Annealing

Simulated annealing was introduced in 1983 by Kirkpatrick, Gelatt, and Vecchi () as a method of function optimization that is particularly well suited to functions that are difficult to evaluate in any continuous manner. It is based on the physical process of annealing, in which the molecules of a material are brought into a crystalline arrangement by gradually cooling the material. Since the crystal is the most ordered configuration of the molecules, it is the one that minimizes the total energy of the system, so the molecules should ultimately come to rest in that configuration. However, there may be other arrangements of the molecules that are stable despite having a higher total energy than the minimum (i.e. local minima of the energy function), and to avoid the system coming to rest in one of these arrangements the scientist must be careful to give the system enough initial energy and to avoid cooling it too quickly. At each stage of the cooling, small temperature fluctuations within the material will create and destroy defects until equilibrium for that temperature is achieved.

Simulated annealing is a procedure for minimization of a function of several variables in which the values of the variables represent configurations of the molecules of the system and the objective function represents its total energy. The computer initially assigns random values to the variables and then gives the system a certain “temperature”, i.e. a tendency of the variables to move randomly. Each move will affect the total energy of the system as measured by the objective function, and the temperature is used to determine the probability that a move that raises the total energy will be accepted. After enough moves have been made to simulate the equilibrium activity of the material at that temperature, the temperature of the system is lowered and the process begins again. Eventually, the temperature is low enough that the variables no longer change significantly and a minimum is achieved. It is important to note, however, that although physical annealing is known in theory to produce the global minimum of the total energy of a system, the efficiency of simulated annealing depends on the function being minimized and it is difficult to guarantee that a given annealing schedule will produce a global minimum rather than a local one.

In our investigation we use an adaptive version of simulated annealing introduced by Lester Ingber (, , , ) as an improvement on his earlier algorithm for Very Fast Simulated Reannealing (, ). ASA takes advantage of structure on the input of the objective function (in our case, the ordering of the $`x_i`$’s described below in Section 4) to decrease the running time of the algorithm and increase the probability that it will find a global minimum of the function. Optimization using VFSR and ASA has been effectively applied in wide variety of situations, including three dimensional image compression (, ), modeling of financial markets (, , , , , , ), dairy farming (, ), neural networks (, , , ), geophysical inversion (), magnetic resonance imaging (), electroencephalography (, ), and combat simulation (, ). An article on the many applications of ASA has appeared in The Wall Street Journal ().

## 4 The Algorithm

The point at which our new method departs from traditional solutions to this problem is in the variable change from the $`w_i`$’s to the $`x_i`$’s in the model (1) above. As we have seen at the end of Section 2, the efficiency of forward and backward regression is not dependent on the ordering of the $`x_i`$’s. In our method, however, we use the following “zipper” algorithm to ensure that if $`i`$ is close to $`j`$, then $`x_i`$ and $`x_j`$ represent similar values of the $`w`$’s: Suppose, for example, that there are only four $`w_i`$’s. Then the $`x_i`$’s would be assigned as follows:

$$\begin{array}{cccc}x_0:=1,\hfill & & & \\ x_1:=w_1,\hfill & x_2:=w_1w_2,\hfill & x_3:=w_1w_3,\hfill & x_4:=w_1w_4,\hfill \\ x_5:=w_4,\hfill & x_6:=w_3w_4,\hfill & x_7:=w_2w_4,\hfill & \\ x_8:=w_2,\hfill & x_9:=w_2w_3,\hfill & & \\ x_{10}:=w_3.\hfill & & & \end{array}$$

(Here we have assumed that there were no categorical characteristics with more than two categories, so that no $`w_i`$ and $`w_j`$ arose from the same characteristic. If they had, then the corresponding cross term $`w_iw_j`$ would be omitted from the assignment above.) By using this algorithm, we attempt to structure the search space in such a way as to make the search for the optimal solution as easy as possible for ASA. To get an idea of why this structuring of the search space is important, consider the difference between trying to find the minimum of a straight line and trying to find the minimum of a series of randomly placed points in the plane. Obviously, the former is much easier (both for a human and a computer). However, the former situation can become the latter very quickly if we allow the points along the $`x`$-axis to shuffle themselves randomly (thereby causing the $`x`$-coordinates of the line to shuffle themselves randomly and removing the ordered structuring of the search space). In the model above, by guaranteeing that $`x_i`$ and $`x_{i+1}`$ always contain exactly one of the $`w_j`$’s in common, we hope that the $`C_p`$ value does not change too radically when transitioning between the two. Of course, we can make no guarantees.

We decide in advance the number $`p`$ of the $`x_i`$’s we want our model to include. This could be any number from 0 to $`n`$, but we found in our simulations with $`n=60`$ that $`p510`$ was usually optimal. This range is based on the outcomes of forward and backward regression, on trial and error with our ASA program, on running time (trials with ten variables took about ten minutes on a Pentium 4 850 MHz processor), and on current practice among researchers in epidemiology.

We then set up dummy variables $`z_i,1ip`$, which may take integer values between 1 and $`n`$ depending on which of the $`x_i`$’s would be included in the model. Our goal is then to find the values of the $`z_i`$’s that minimize the objective function defined by the $`C_p`$ statistic for a given model. This is where we use simulated annealing.

The publicly available C-code for ASA () allows the user to define how to evaluate the objective function to be minimized. In our case, evaluation on a particular set of values of the $`z_i`$’s required building a model with the corresponding $`x_i`$’s included, finding the coefficients $`\beta _i`$ by calculating a least squares solution to $`\text{X}\stackrel{}{\beta }=\stackrel{}{y}`$, and calculating the $`C_p`$ statistic for that model. We modified the public code accordingly and ran multiple simulations with various sizes of data sets.

## 5 The Results

We ran several simulations on computer-generated data. In all cases, we used our ASA method and also ran forward and backward regression on the data as a control. We compared the three methods according to the final $`C_p`$ statistic produced and the running time required to achieve it as measured by the number of function evaluations involved. Here one function evaluation is defined to be the calculation of the $`C_p`$ statistic for a given set of variables as described at the end of Section 4, since forward and backward regression use the same procedure as they search for variables to be added to or deleted from the model. It is true that not all function evaluations will take equal amounts of processor time, since an evaluation with more variables will involve finding the pseudoinverse of a larger matrix; however, we believe, in keeping with standard practice in computer science, that such a measure is still meaningful enough to merit study.

We show below the results of our simulations on computer-generated data. Using the characteristics described above in Section 2, there were 60 possible variables; we set our ASA program to search for the best seven variables. We ran five simulations each on populations of size 100,000, 500,000, and 1,000,000 with the following results.

Simulated populations of 100,000 (approximately 45 diseased, 45 healthy studied).

|  | Forward Regression | Backward Regression | ASA |
| --- | --- | --- | --- |
| $`\begin{array}{c}C_p\text{ statistic}\hfill \end{array}`$ | $`\begin{array}{c}46.387\hfill \\ 45.831\hfill \\ 54.046\hfill \\ 56.199\hfill \\ 54.936\hfill \end{array}`$ | $`\begin{array}{c}29.574\hfill \\ 43.318\hfill \\ 50.251\hfill \\ 56.199\hfill \\ 49.022\hfill \end{array}`$ | $`\begin{array}{c}46.734\hfill \\ 48.968\hfill \\ 53.009\hfill \\ 56.199\hfill \\ 54.724\hfill \end{array}`$ |
| $`\begin{array}{c}\text{Function\hspace{0.17em} evaluations}\hfill \end{array}`$ | $`\begin{array}{c}177\hfill \\ 119\hfill \\ 177\hfill \\ 119\hfill \\ 177\hfill \end{array}`$ | $`\begin{array}{c}1725\hfill \\ 1802\hfill \\ 1824\hfill \\ 1830\hfill \\ 1820\hfill \end{array}`$ | $`\begin{array}{c}1341\hfill \\ 1677\hfill \\ 1587\hfill \\ 2827\hfill \\ 199\hfill \end{array}`$ |

Simulated populations of 500,000 (approximately 225 diseased, 225 healthy studied).

|  | Forward Regression | Backward Regression | ASA |
| --- | --- | --- | --- |
| $`\begin{array}{c}C_p\text{ statistic}\hfill \end{array}`$ | $`\begin{array}{c}47.511\hfill \\ 43.141\hfill \\ 50.803\hfill \\ 36.406\hfill \\ 45.648\hfill \end{array}`$ | $`\begin{array}{c}41.831\hfill \\ 37.691\hfill \\ 41.506\hfill \\ 35.767\hfill \\ 41.693\hfill \end{array}`$ | $`\begin{array}{c}47.510\hfill \\ 43.152\hfill \\ 51.160\hfill \\ 41.127\hfill \\ 46.063\hfill \end{array}`$ |
| $`\begin{array}{c}\text{Function\hspace{0.17em} evaluations}\hfill \end{array}`$ | $`\begin{array}{c}119\hfill \\ 399\hfill \\ 290\hfill \\ 119\hfill \\ 452\hfill \end{array}`$ | $`\begin{array}{c}1815\hfill \\ 1785\hfill \\ 1764\hfill \\ 1739\hfill \\ 1752\hfill \end{array}`$ | $`\begin{array}{c}2442\hfill \\ 3336\hfill \\ 1305\hfill \\ 3780\hfill \\ 2994\hfill \end{array}`$ |

Simulated populations of 1,000,000 (approximately 450 diseased, 450 healthy studied).

|  | Forward Regression | Backward Regression | ASA |
| --- | --- | --- | --- |
| $`\begin{array}{c}C_p\text{ statistic}\hfill \end{array}`$ | $`\begin{array}{c}47.162\hfill \\ 44.149\hfill \\ 40.755\hfill \\ 16.518\hfill \\ 43.012\hfill \end{array}`$ | $`\begin{array}{c}39.287\hfill \\ 41.134\hfill \\ 21.738\hfill \\ 22.845\hfill \\ 40.003\hfill \end{array}`$ | $`\begin{array}{c}48.347\hfill \\ 44.148\hfill \\ 40.743\hfill \\ 46.679\hfill \\ 38.235\hfill \end{array}`$ |
| $`\begin{array}{c}\text{Function\hspace{0.17em} evaluations}\hfill \end{array}`$ | $`\begin{array}{c}234\hfill \\ 345\hfill \\ 234\hfill \\ 177\hfill \\ 504\hfill \end{array}`$ | $`\begin{array}{c}1764\hfill \\ 1764\hfill \\ 1659\hfill \\ 1620\hfill \\ 1752\hfill \end{array}`$ | $`\begin{array}{c}461\hfill \\ 329\hfill \\ 678\hfill \\ 690\hfill \\ 1913\hfill \end{array}`$ |

Asymptotically, the number of function evaluations necessary to compute forward regression is easily seen to be $`O(nk)`$ where $`n`$ represent the total number of variables under consideration and $`k`$ represents the final number of variables chosen in the model. It is a notoriously difficult problem of theoretical computer science to analyze the computational complexity of ASA . We can only present our results and leave it to the interested reader to draw his own conclusions.

## 6 Advantages and Limitations

We found that running ASA in most cases produces a slightly lower $`C_p`$ statistic especially on the trials with large data sets. The cases in which ASA did not produce a lower $`C_p`$ statistic were ones in which the optimal solution required significantly more or fewer than the seven variables we programmed ASA to search for; in these cases the $`C_p`$ statistic produced by ASA was still close to that of the other methods, and when we reconfigured ASA to search for an appropriate number of variables it produced a lower $`C_p`$ statistic. It appears that for small data sets with few parameters, the traditional methods of modeling are slightly preferable to ASA. The benefits of ASA seem to be more pronounced for larger data sets.

Another advantage of ASA is that its results seem to be more predictable than those of the other methods. Although in many tests the results of the three methods were comparable, in several instances forward regression and backward regression produced $`C_p`$ statistics that were considerably higher than the lowest obtained $`C_p`$ statistic; the statistics produced by ASA were never far from the best statistics produced by the three methods.

One disadvantage of our method is that it requires the user to determine the number of variables in the model beforehand. Simulations run under the same conditions can have optimal solutions (as determined by forward and backward regression) requiring quite different numbers of variables; in our simulations with 60 possible nonconstant variables we found some for which the $`C_p`$ statistic was optimized with just one variable and some requiring as many as 21.

In practice this should not be a serious concern for three reasons. One is that a researcher in epidemiology will usually have a rough idea of the desired number of variables in advance as a compromise between the needs of the problem and the availability of computing power. (The assumption above that the model is linear already requires that the researcher know in advance something about the desired model.) A second reason is that as part of the random component of the program it automatically checks the possibility of deleting variables from the model, so if the minimum involves fewer variables than the user selects, the program will still find it. Finally, this property may even be considered an advantage, since our method allows the user to control the sophistication of the model (i.e. the number of variables it will use) in advance, whereas with forward and backward regression the models produced may have widely varying numbers of variables.

Another disadvantage of our method is that it is not guaranteed to find the global minimum of the function we are trying to optimize. Neither, however, do the traditional methods of regression, and indeed, such a goal would be impractical because of the number of parameters in the problem.

## 7 Further Research

Though the results of the computer simulations are suggestive, we would like to see more done in the area of practical epidemiology, perhaps a large-scale analysis of publicly available data sets. Such an undertaking is non-trivial, however, even given this previous research. Most publicly available data sets are enormous in comparison with the relatively small ones used in this study (several hundreds versus over fifteen thousand). Also, the number of possible choices for variables is immense. The investment of computer time would be substantial even for a modestly ambitious study. An even greater hurdle seems to be the numerical stability of the algorithm itself. For data sets that are extremely large or that have a very wide range of data values, finding the least squares solution (which involves inverting a matrix) is inherently a numerically unstable undertaking. We have attempted programming the bulk of the algorithm in both C and Mathematica. In C, the numerical stability of the unusually large number of computations seems to be an issue; whereas, in Mathematica, even though the kernel controls the numerical precision well, the time necessary for even a relatively small computation is enormous.


text:
# 1

## 1

In recent years, the effective field theory (EFT) method has become the main tool to deal with various low energy processes and strong interactions in the nonperturbative regime. Beginning from Weinberg’s seminal works, this method has been extensively and successfully applied to the low energy nucleon systems. As an EFT parametrizes the high energy details in a simple way, there appear severe UV divergences (or ill-definedness) and more undetermined constants that must be fixed somehow. Owing to the nonperturbative nature, the EFT for nuclear forces also becomes a theoretical laboratory for studying the nonperturbative renormalization and a variety of renormalization (and power counting) schemes have been proposed*** We apologize in advance for those contributions we have not been informed of by now.. However, a complete consensus and understanding on this important issue has yet not been reached. Though the EFT philosophy is very natural, but its implementation requires discretion (in parametrizing the short distance effects). In particular, as will be shown below, the nonperturbative context impedes the implementation of conventional subtraction renormalization, and some regularization/renormalization (R/R) prescription might even fail the EFT method or physical predictions, unlike in the perturbative case.

In this short report, we shall revisit the problem without specifying the R/R scheme (we only need that the renormalization is done) to elucidate the key features for the implementation of R/R in nonperturbative regime, and various proposals could be compared and discussed. Then the key principles that should be followed are suggested as the first steps in such understandings.

## 2

The physical object we are going to investigate is the transition matrix $`T`$ for nucleon-nucleon scattering processes at low energies,

$$N(𝐩)+N(𝐩)N(𝐩^{})+N(𝐩^{}),$$
(1)
which satisfies the Lippmann-Schwinger (LS) equation in partial wave formalism,

$`T_{ll^{}}(p^{},p;E)=V_{ll^{}}(p^{},p)+{\displaystyle \underset{l^{\prime \prime }}{}}{\displaystyle \frac{kdk^2}{(2\pi )^2}V_{ll^{\prime \prime }}(p^{},k)G_0(k;E^+)T_{l^{\prime \prime }l^{}}(k,p;E)},`$
$`G_0(k;E^+){\displaystyle \frac{1}{E^+k^2/(2\mu )}},E^+E+iϵ,`$ (2)
with $`E`$ and $`\mu `$ being respectively the center energy and the reduced mass, $`𝐩`$ ($`𝐩^{}`$) being the momentum vector for the incoming (outgoing) nucleon, and $`p^{}=|𝐩^{}|`$, $`p=|𝐩|`$. The potential $`V(p^{},p)`$ can be systematically constructed from the $`\chi `$PT. We remind that the constructed potential is understood to be finite first, as ’tree’ vertices in the usual field theory terminology.

Now let us transform the above LS equation into the following compact form as a nonperturbative parametrization of $`T`$-matrix, we temporarily drop all the subscripts that labelling the status of the finiteness and divergence:

$`T^1=V^1𝒢,`$ (3)
where $`𝒢V^1\{\frac{kdk^2}{(2\pi )^2}VG_0T\}T^1`$ is obviously a nonperturbative quantity. Here $`V`$, $`T`$ and therefore $`𝒢`$ are momenta dependent matrix in angular momentum space. In field theory language, $`𝒢`$ contains all the ’loop’ processes generated by the 4-nucleon vertex $`V`$ joined by nucleon (and antinucleon) lines. This nonperturbative parametrization clearly separates $`V`$ (Born amplitude) from $`𝒢`$ (all the iterations or rescatterings), and on-shell ($`p^{}=p=\sqrt{2\mu E}`$) unitarity (nonperturbative relation) of $`T`$ is automatically satisfied here thanks to the on-shell relation between the $`K`$-matrix and $`T`$-matrix

$`T_{os}^1=K_{os}^1+i\mu {\displaystyle \frac{\sqrt{2\mu E}}{2\pi }},`$ (4)
Using Eq.( 3) we have

$`𝒢_{os}=V_{os}^1K_{os}^1i\mu {\displaystyle \frac{\sqrt{2\mu E}}{2\pi }}.`$ (5)
Any approximation to the quantity $`𝒢`$ is leads to a nonperturbative scheme for $`T`$.

Since the on-shell $`T`$-matrix must be physical, then for general potential functions this nonperturbative parametrization leads to the following obvious but important observations in turn: (I) it is impossible to renormalize $`𝒢`$ through conventional perturbative counter terms in $`V`$ (we will demonstrate this point shortly), i.e., they have to be separately renormalized; (II) the perturbative pattern of the scheme dependence for transition matrix or physical observables breaks down in the nonperturbative contexts like Eq.( 3); (III) the R/R for $`V`$ and $`𝒢`$ must be so performed that the on-shell $`T`$-matrix satisfies physical boundary conditions (bound state and/or resonance poles, scattering lengths or phase shifts, etc.), i.e, $`V`$ and $`𝒢`$ are ’locked’ with each other; (IV) due to this ’locking’, the validity of the EFT power counting in constructing $`V`$ also depends on the R/R prescription for $`𝒢`$, i.e., the power counting of $`V`$ could be vitiated by the R/R prescription for $`𝒢`$ which is performed afterwards. This explains why there are a lot of debates on the consistent power counting schemes for constructing $`V`$. Point (III) is in fact what most authors actually did or is implicit in their treatments. we will return to these points later.

To be more concrete or to demonstrate the first point (observation (I)), let us consider the case without mixing ($`l=l^{}`$) for simplicity. Then the matrices become diagonal and only energy dependent due to on-shell condition, we have the following simple form of nonperturbative parametrization (from now on we will drop the subscript ”os” for on-shell)We can also arrive at it simply through dividing both sides of Eq.( 2) by $`V_l(E)`$ and $`T_l(E)`$ and rearranging the terms.,

$`{\displaystyle \frac{1}{T_l(E)}}={\displaystyle \frac{1}{V_l(E)}}𝒢_l(E),`$ (6)
where $`𝒢_l(E)=\frac{({\scriptscriptstyle {\scriptscriptstyle \frac{kdk^2}{(2\pi )^2}}VG_0T})_l(E)}{V_l(E)T_l(E)}`$. In perturbative formulation, all the divergences are removed order by order before the loop diagrams are summed up. Here, in Eq.( 6) (or Eq.( 3)), the renormalization are desired to be performed on the compact nonperturbative parametrization. After $`V`$ is calculated and renormalized within $`\chi `$PT, the task is (a) to remove the divergence or ill-definedness in $`𝒢_l(E)`$ and (b) to make sure that the renormalized quantities ($`V_l^{(R)}(E)`$ and $`𝒢_l^{(R)}(E)`$) lead to a physical $`T`$-matrix. The most important one is (b), namely the $`T`$-matrix obtained must possess reasonable or desirable physical behaviors after removing the divergences. Technically, this amounts to a stringent criterion for the R/R prescription in use: the functional form of $`𝒢_l(E)`$ could not be altered, only the divergent parameters (coefficients) get ’replaced’ by finite ones.

Now let us suppose $`𝒢_l`$ could be renormalized by introducing counter terms. Within this parametrization of Eq.( 6), the potential $`V`$ appears as the only candidate to bear counter terms, then from Eq.( 6) we would have

$`𝒢_l^{(R)}(E)=𝒢_l^{(B)}(E)+{\displaystyle \frac{1}{V_l^{(R)}(E)}}{\displaystyle \frac{1}{V_l^{(R)}(E)+\delta V_l(E)}},`$ (7)
where the superscripts $`(R)`$ and $`(B)`$ refer to ’renormalized’ and ’bare’, $`\delta V_l`$ denotes the additive counter terms and $`V_l^{(B)}(E)V_l^{(R)}(E)+\delta V_l(E)`$.

From the definition below Eq.( 3), $`𝒢`$ must take the following form which is no less complicated than a nonpolynomial function in $`p(=\sqrt{2\mu E})`$,

$$𝒢_l^{(B)}(E)=\frac{\frac{d^3k}{(2\pi )^3}V_l(p,k)G_0(k;E^+)T_l(k,p;E)}{V_l(E)T_l(E)}=\frac{N_l^{(B)}(p)}{D_l^{(B)}(p)}\frac{\mu }{2\pi }ip,$$
(8)
with $`N_l^{(B)}(p),D_l^{(B)}(p)`$ being at least two nontrivial divergent or ill-defined polynomials. With this parametrization, Eq.( 7) becomes

$$𝒢_l^{(R)}(E)=\frac{N_l^{(B)}(p)V_l^{(R)}(E)V_l^{(B)}(E)+D_l^{(B)}(p)\delta V_l(E)}{D_l^{(B)}(p)V_l^{(R)}(E)V_l^{(B)}(E)}\frac{\mu }{2\pi }ip=\frac{N_l^{(R)}(p)}{D_l^{(R)}(p)}\frac{\mu }{2\pi }ip.$$
(9)
Note that in the divergent or ill-defined fraction $`\frac{N_l^{(B)}(p)V_l^{(R)}(E)V_l^{(B)}(E)+D_l^{(B)}(p)\delta V_l(E)}{D_l^{(B)}(p)V_l^{(R)}(E)V_l^{(B)}(E)}`$ only $`V_l^{(R)}(E)`$ is finite. Then except for the simplest form of $`V_l(E)`$ (= $`C_0`$, see, e.g., Ref.), it is impossible to obtain a nontrivial and finite fraction out of the this divergent fraction with whatever counter terms $`\delta V_l(E)`$, i.e., it is impossible to remove the divergences in the numerator and the denominator at the same time with the same counter terms as in Eq.( 9). Examples for nonperturbative divergent fractions could be found in Ref.. Even by chance $`𝒢`$ is finite, then its functional form in terms of physical parameters ($`E,p`$) must have been altered (often oversimplified) after letting the cutoffs go infinite. For channels with mixing, with $`V_l`$ and $`𝒢_l`$ becoming $`2\times 2`$ matrices, all the preceding deductions still apply. Thus within nonperturbative regime the counter term (via potential) renormalization of $`T`$-matrix fails, that is, $`𝒢`$ must be separately renormalized–the first observation given above is demonstrated as promised. The other observations follow easily.

What we have just shown does not mean that the counter term fails at any rate. The above arguments showed that the counter term could not successfully renormalize the nonperturbative $`𝒢`$ or $`T`$-matrix except they are introduced before the corresponding infinite perturbative series are summed up or before all parts of the nonperturbative objects are put together and ’installed’, say, before the corresponding Schrödinger equation is solved. It is well known that the Schrödinger equation approach or similar is intrinsically nonperturbative, therefore the counter terms introduced there (as parts of the potential operators) ARE nonperturbative in the sense that they enter before the $`T`$-matrix or phase shifts are finally derived, or equivalently, they must ’enter’ before the perturbative series are summed up. So there is no contradiction between our conclusion here and that in Ref.. In fact they are in complete harmony and our conclusion above turns out to be also supporting the Schrödinger equation approach for renormalizing nucleon-nucleon scattering. Thus our discussions above are helpful to clarify the difference between perturbative and nonperturbative renormalization or to remove some mists or controversies around the renormalization of EFT for nucleon-nucleon scattering.

To make our language more concise, we would temporarily call the counter terms ’endogenous’ to an object if they must enter before all the necessary calculations on the interested object are done, and ’exogenous’ otherwise.

## 3

Therefore the renormalization in nonperturbative regime have to be implemented either with ’endogenous’ counter terms (e.g.,) or otherwise (e.g., like in Refs.). To proceed we recall the conceptual foundation for EFT method: there must be a complete theory underlying all the low energy EFTs. Suppose we could compute the $`T`$-matrix in the underlying theory, then in the low energy limit, we should obtain a finite matrix element parametrized in terms of EFT coupling constants and certain finite constants arising from the limit. It is these constants that we are after. In the ill-defined EFT frameworks, a R/R or subtraction prescription is employed to ’retrieve’ the finite constants. As we have seen the failure of the counter term subtraction at the potential or vertex level in nonperturbative parametrization, the renormalization must be otherwise implemented. A prompt choice is the subtraction at the integral level, and the counter term formalism is dismissed from the underlying theory point of view. In this sense, the formal consistency issue of Weinberg power countings is naturally dissolved without resorting to other means, and the real concern is the renormalization instead of power counting rules. This resolution is in fact just a paraphrase of the fourth observation in section two (point (IV)): a consistent power counting could be vitiated by the afterward renormalization procedures such as those using exogenous counter terms.

If the ill-definedness only originates from the local part of the potential, the integral level subtraction is OK as the divergent integrals can be cleanly isolatedFor more general parametrization of the ill-defined integrals, see.. For the more interesting cases with general pion-exchange contributions, such clean isolation seems impossible, which leads some authors to perturbative approach (KSW), but the convergence is shown to be slow and the KSW power counting also bears some problems. Due to such difficulties, the regularization schemes with finite cutoffs or parameters have also been adopted to keep the investigations nonperturbative§§§All the finite cutoff regularization schemes could be seen as effectively employing endogenous counter terms implicitly: the UV region contributions are subtracted before integration., which have to be numerical. In such approaches certain sort of ’fitting to data’ is incorporated in one way or the other, this is what we stressed in section two. There recently appears a new approach that constructs the nonperturbative $`T`$-matrix from the KSW approach, where again certain sort of fitting is crucial there.

Then it is clear that the main obstacle is the severe regularization dependence in nonperturbative regime. If we could parametrize all the ill-definedness in ambiguous but finite expressions, then the obstacle would disappear and the only task is to fix the ambiguities. Otherwise, any theoretical judgement just based on one regularization scheme is unwarranted or ungrounded unless it is justified in various schemes . The importance of regularization in nonperturbative regime could also be seen in the following way: For any R/R prescription in nonperturbative regimes to be able to ’reproduce’ such ’true’ $`T`$-matrix obtained from underlying theory, the regularization in use should at least reproduce the same nonperturbative functional expression as that given in underlying theory, with certain constants being adjustable, otherwise the ’true’ functional expression could not be approached since exogenous counter terms, as we have shown above, could not be effected in nonperturbative parametrization. Alternatively, the foregoing arguments could be interpreted as a call for a devise of a new nonperturbative framework where endogenous counter terms could be naturally incorporated.

Our parametrization also points towards a new technical direction for the nonperturbative determination of the $`T`$-matrix: to focus on the calculation of $`𝒢_l(E)`$. The relevance of R/R in $`𝒢`$ can be made clearer in the low energy expansion of $`𝒢`$ (Padé or Taylor) since it is necessarily a function of the center energy ($`E`$) or on-shell momentum ($`p`$):

$$Re(𝒢_l(p))|_{\text{Padé}}=\frac{n_{l;0}+n_{l;1}p^2+\mathrm{}}{d_{l;0}+d_{l;1}p^2+\mathrm{}}|_{\text{Taylor}}=g_{l;0}+g_{l;1}p^2+g_{l;2}p^4+\mathrm{}.$$
(10)
Here the coefficients $`g_{l;n}`$ will inevitably be R/R dependent and must be determined together with the parameters in the potential from physical boundary conditions, say the empirical phase shifts in the low energy ends. Investigations along this line are in progress and the primary the results are interesting. In fact, there is several virtues in such analysis: Firstly, one needs not to carry out the concrete the renormalization of $`𝒢`$; Secondly, by comparing with low energy region of the empirical data of $`NN`$ scattering (say, $`E50\text{Mev}`$) one could in principle easily obtain the optimal values for the semi-phenomenological parameters ($`g_{l;n}`$) which in turn could yield a pretty good prediction of higher energy behaviors though the parametrization formula of Eq.( 3) is strikingly simple; Thirdly, with such simple tools one could test whether EFT systematically works for $`NN`$ scattering; Fourthly, one could also use it to estimate if the potential constructed up to a certain order is sufficient for various purposes without any loop calculations and renormalization operations. Moreover, due to its simplicity in principle, one might also find other uses like the estimation of the coupling constants of the local terms in constructed potential. This list could be further enlarged. We will visit them in the near future. In Fig. 1 we serve a primary example for such analysis employing the potential given by EGM: the phase shifts of $`{}_{}{}^{1}D_{2}^{}`$ channel predicted only with a single parameter $`Re(𝒢_{l=2})g_{l=2;0}`$ with potentials input at leading order (Lo), next-to-leading order (Nlo) and next-to-next-to-leading order (Nnlo), respectively. It is clear that the prediction of the phase shift for the range $`50\text{Mev}E200\text{Mev}`$ with only one optimized or fitted parameter of $`Re(𝒢_{l=2})`$, $`g_{l=2;0}`$, improves with chiral order for the construction of the potential (Lo, Nlo, and Nnlo). For most channels this seems true. Of course, there could well be exceptions, which we interpret it as insufficiency of the chiral expansion in the corresponding channel(s).

In summary, we proposed a simple and novel parametrization for understanding the nonperturbative renormalization of the $`T`$-matrix for nucleon-nucleon scattering. The distinctive feature of the nonperturbative renormalization–the failure of ’exogenous’ counter terms or the need of ’endogenous’ ones is clarified together with some related issues like consistency of W power counting. By the way we could derive other uses from our simple parametrization of the $`T`$-matrix that might be of some academic values. We hope our investigation might be useful for a number of important issues related to nucleon interactions, and also to other nonperturbative problems.

## Acknowledgement

I wish to thank Dr. E. Ruiz Arriola and Dr. J. Gegelia for helpful communications on the topic. I am especially grateful to an anonymous referee for his very helpful and encouraging report that leads to the present form of this manuscript. This project is supported in part by the National Natural Science Foundation under Grant No. 10502004.


text:
# Zipf Law in Firms Bankruptcy

## 1 Introduction

It is a challenging task to model aggregate dynamics of firms in econophysics. As such dynamics for firms with considerable amount of money flow is concerned, one can now have evidence of several phenomenological facts about the dynamics of firms , which are not only interesting in physics, but also have quite significance in macroeconomics .

Among those facts are size distribution and growth of firms. Distribution of firm size such as sales, number of employees, total assets follows a power-law, or Pareto-Zipf law . As for the growth of firms, the distribution of the logarithm of growth rates has an exponential form and that the fluctuations in the growth rates scale with the firm size . While these observation concerns about stock of a firm, no less important is flow such as profits, whose properties of distribution and growth rates have recently observed in empirical works .

While these phenomenological facts have been found for living firms, little is known about the death in aggregate dynamics of a large number of firms. That is exactly what the present work concerns about. It is important to uncover phenomenology of how firms die in order to understand the firms growth dynamics.

In section 2, we make observation of Japanese bankruptcy by using an exhaustive list of bankrupted firms in a year. We show phenomenological findings about the distribution of total liabilities or debts when bankrupted, the life-time of bankrupted firms, and correlation between the firm’s debt and size. These are the main results of this paper. In section 3, we describe a firm’s financial activity in terms of balance-sheet dynamics. A model for debtor of firms and creditor of a bank, recently proposed by economists , after briefly recapitulated, is shown to be able to reproduce our findings about bankruptcy. Summary and discussion is given in section 4. In appendix, we give description about the nature of bankruptcy dataset employed in this paper.

## 2 Firms Bankruptcy in Japan

Our data, provided by a domestic data-bank company (Tokyo Shoko Research, Ltd.), is an exhaustive list of Japanese bankrupted firms in 1997. The data is exhaustive in the sense that any bankrupted firm with total amount of debt exceeding 10 million yen (roughly 80,000 euros/dollars) is listed in it. The number of such firms in the year is 16,526. The data contains information of debt, sales, business sector, number of employees, primary bank when bankrupted, dates of establishment and failure, and stockholders fund when established. Two remarks immediately follow.

First, bankruptcy or business failure is not a legal term, but should be understood as a critical financial insolvency of a debtor. Typical case and classification of different types are given in appendix. What is important for our purpose here, however, is that most of the cases in the dataset, whatever the types are, is basically caused by financial insolvency state of the firm. Indeed, out of the total 16,526 cases, the most frequent one is suspension of bank transactions (13,850). This liquidity problem (see appendix for details) is considered as the main cause for other legal and private procedures of bankruptcy.

Second, debt refers to the total sum of liabilities in the credit (“right-hand”) side of balance-sheet. This quantity, when the firm goes bankrupted, is relatively easy to measure in comparison with other balance-sheet quantities such as total asset and equity111Debt here is not net liabilities defined by the total liabilities minus the amount of assets not hypothecated, which is difficult to know.. Since each item in debt has a creditor outside of the firm, so the quantity is not basically hidden inside. Actual measurement of the amount of debt is based either on legal documents or on investigation of current/non-current liabilities in balance-sheet plus the amount of bills discounted and transfers by endorsement, according to legal and private procedure for bankruptcy, respectively.

We used the entire dataset in what follows.

### 2.1 Distribution of debt

Fig. 1 shows the cumulative probability distribution of debt when the firm was bankrupted. For large debts exceeding $`10^8`$ yen, we can observe a power-law distribution over three orders of magnitude or even more. The probability that a given bankrupted firm has debt equal to or greater than $`x`$, denoted by $`P_>(x)`$, obeys

$$P_>(x)x^\mu ,$$
(1)
with a constant $`\mu `$, the so-called Pareto index for bankruptcy. The index was estimated as $`\mu =0.911\pm 0.008`$ (by least-square-fit for samples equally spaced in logarithm of rank in the power-law region; error is 90% level). This phenomenological finding is what we call Zipf-law for debt of firms bankrupted. It was pointed out in  that debt may obey Zipf-law, but unfortunately with a quite limited number of data (corresponding to the tail more than $`10^{11}`$ yen of debt).

### 2.2 Life-time of bankrupted firms

No less interesting than how large a firm can fail is the problem of how long a failed firm lived. Fig. 2 gives the cumulative distribution of life-time $`\tau `$, $`P_>(\tau )`$, for the same set of bankrupted firms. With this plot being semi-log, it is observed that the entire distribution follows an exponential distribution. We also see that the distribution has a sudden drop and kink in the shape.

In order to quantify the shape of $`P_>(\tau )`$, let us introduce a function defined by

$$h(t)=\frac{\frac{dP_>(t)}{dt}}{P_>(t)}.$$
(2)
It is easily seen from the right-hand side expression of eq. (2) that for a sufficiently small time-interval $`\mathrm{\Delta }t`$, $`h(t)\mathrm{\Delta }t`$ is the probability that the life-time $`\tau `$ satisfies $`t<\tau <t+\mathrm{\Delta }t`$ under the condition that $`\tau t`$. This function has a similar notion as hazard rate function in probability (see  for example). If $`h(t)`$ is constant $`h_0`$, then $`P_>(t)=\mathrm{exp}(h_0t)`$. $`h(t)`$ has the physical dimension of inverse of time.

From $`P_>(t)`$ in Fig. 2, we estimate the derivative in eq. (2) by doing moving-average of stepwise slopes over one year. The resulting $`h(t)`$ is given in Fig. 3 (a). $`h(t)`$ is nearly constant over 30 years, but takes large values for about 10 years right after the World War II. In general, the number of firm $`t`$ years old can be estimated by multiplying the number of firms born $`t`$ years ago by the probability of survival of a firm to an age of $`t`$ years. Therefore, if the process of survival is homogeneous in time, the time-scale in $`h(t)`$ comes from entry rate of new firms. Fig. 3 (b) is the historical record of entry of firms . One can observe that the 10 years period of high $`h(t)`$ corresponds to the period of extremely high birth-rate right after the world war II for about 10 years. This was brought about by the dissolution of zaibatsu (corporate alliances) monopoly under the government by U.S. General Headquarters (from 1945 to 1952). One can also note that $`h(t)`$ has the same order of magnitude as the entry rate of firms by converting time from days to years all over the history.

### 2.3 Correlation between firm size and debt when bankrupted

Since it has been known222But see also another paradigm . that firm size follows Zipf law , let us take a look at relation between firm size and firm debt right before bankruptcy. Fig. 4 shows a clear correlation between them. Here the firm size is measured by sales right before bankruptcy. The number of employees, as a proxy of firm size, gives similar result as for sales. The bigger a firm is, the larger its debt is when failed.

To quantify this relation, denoting sales of a bankrupted firm by $`s`$, we use the ratio $`A=x/s`$ where debt $`x`$ and sales $`s`$ are given in the same units of money. Introduce the logarithm of $`A`$ as $`a=\mathrm{log}_{10}(A)`$, and we examined the probability density of $`a`$, $`P(a)`$, in Fig. 5 (filled dots) by using all the data $`x>10^8`$ (yen), which corresponds to the Pareto-Zipf range in eq. (1). It is found that the probability density function $`P(a)`$ is exponential distribution with skewness, i.e., with different slopes for $`a>0`$ and $`a<0`$. The side of $`a>0`$ is fatter than the other side because a positive $`a`$ means that the firm is nearer to the state of insolvency.

Furthermore, we show that the ratio does not depend statistically on which range of debt is being observed. Preparing bins as $`x[10^{8+0.4(n1)},10^{8+0.4n}]`$ yen ($`n=1,\mathrm{},4`$), we can examine the conditional probability density, $`P(a|x)`$, for each bin of $`x`$. Fig. 5 shows that statistical dependence of the ratio on the value of debt is very weak. This means that the relation between firm size and debt is multiplicative with multiplying factor $`A`$ drawn from the distribution given in Fig. 5. If the statistical independence between $`A`$ and $`s`$ holds over the entire range of the variables, the probability density for $`\mathrm{log}x`$ is equal to the convolution of the one for $`a=\mathrm{log}A`$ and the one for $`\mathrm{log}s`$. If $`s`$ has a power-law distribution with a Pareto index, then $`x`$ will have a power-law with the same Pareto index. Though being a non-rigorous argument because of the breakdown of power-law under some threshold and of the statistical independence, this give us a basic idea that the Zipf law for firm size and that for firm debt when bankrupted are closely related to one another.

## 3 Balance-sheet dynamics of firm

### 3.1 What is a firm?

In order to understand how firms die, let us consider the activity of a firm from the viewpoint of money flow.

A firm makes activity by producing goods or services in anticipation to be demanded by consumers, and to yield profit. For the firm to realize the production and selling of goods and services, it must usually invest more money than what it actually possesses with a fund by the firm’s owners or stockholders. The investment is done not only for cost of sales, expenses of selling, etc. but also for additional facility, production line, employment etc.

A firm’s activity has three facets in general:

1. finance
borrowing money from banks, market investors
2. investment
spending money in anticipation of future return of profits
3. collecting
earning profits from the sales of good and services

In other words, these three facets imply (1) the presence of creditor (borrower), (2) that investment can be a source of risk, because (3) return of profits has uncertainty in realization.

The above description can be formulated in terms of stock and flow of money. Flow is inflow of money in a period of time, and stock is the accumulation of flow over time. Stock of a firm is its balance sheet which consists of asset, debt and equity (see Fig. 6). Here asset refers to all the assets possessed by the firm which are expected to yield profits. The asset can be classified into equity and debt according to how it is financed, either by stockholders or by loan from banks and market, respectively.

Temporal change of balance-sheet is determined by flow. One source of flow is profit which is basically given by sales minus cost of sales, finance, and others. The other flow is investment which is determined by the firm’s decision based on anticipation of profits. Since the investment is covered by additional debt, it also depends on creditor’s activity.

Balance sheet dynamics is the temporal change of asset and debt. Debt is necessary for making profits, but is a source of risk since it can diminish equity and can dominate the balance sheet. Financial fragility can be measured either by flow variable such as insolvency (the ratio between debt commitments to profits) or by stock variable such as liquidity (the ratio between debt and capital). Persistent state of insolvency results in liquidity problem. Actually the most cases of bankruptcies (see Appendix) are due to the problem of liquidity.

### 3.2 Aggregate dynamics of creditor and debtors

Recently a model for the dynamics of balance sheets of a bank and firms are proposed in , which is based on their previous model for aggregate dynamics of firms . This model of creditor-debtor dynamics has two elements; competition between firms by their sizes, and debtor-creditor relation between a bank and each firm.

In this section, let us briefly recapitulate the model (see  for details) and show by simulation that our phenomenological findings about bankruptcy can be reproduced in the model.

A firm $`i`$ has, at time $`t1`$, asset $`K_{t1}^i`$, debt $`L_{t1}^i`$ and equity $`A_{t1}^i`$ which satisfy the identity:

$$K_{t1}^i=L_{t1}^i+A_{t1}^i.$$
(3)
$`i=1,\mathrm{},n`$ where $`n`$ is the number of firms. In the period of time from time $`t1`$ to $`t`$, the firm makes the following activity with $`K_t^i`$ to be determined below. It produces an output $`\varphi K_t^i`$ ($`0<\varphi <1`$), which is to yield revenue with uncertainty. At the same time, the firms has financial cost, comprised of dividend for equity and interest rate for debt. The profit in the current period is

$$\pi _t^i=u_t^i\varphi K_t^ir_t^iK_t^i,$$
(4)
where $`r_t^i`$ is interest-rate for the firm $`i`$. Here $`u_t^i`$ is a stochastic variable representing the uncertainty in making profit. In the second term, the ratio of dividend to equity is assumed to be equal to the interest rate.

Supposing $`r_t^i`$ is a parameter, the firm plans the amount of investment during the current period. Since the profit $`\pi _t^i`$ is proportional to the total asset $`K_t^i`$, the firm could increase its size by investment. But the bigger the size is, the larger the risk of insolvency or bankruptcy is, in the presence of uncertainty for profit. As a simplest form, the objective function for determining the investment is assumed to have a quadratic term of $`K_t^i`$. Namely, the maximization of the expectation:

$$E[\pi _t^ic(K_t^i)^2]$$
(5)
yields an optimal value of $`K_t^i`$. Then the investment in the current term is

$$I_t^i=K_t^iK_{t1}^i$$
(6)
for the realization of which the firm demand loan of the amount

$$L_t^{i(d)}=L_{t1}^i\pi _{t1}^i+I_t^i.$$
(7)
In this model the demand is assumed to be always fulfilled by bank with an appropriate interest-rate. Therefore balance sheet updates by $`K_t^i=K_t^i`$, $`L_t^i=L_t^{i(d)}`$, and $`A_t^i=K_t^iL_t^i`$, once the interest rate $`r_t^i`$ is determined as shortly given.

Bank is monopolistic and responsible for all the loans to firms. Its balance sheet is composed of loan $`L_{t1}`$ as asset, deposit $`D_{t1}`$ as debt, and equity $`E_{t1}`$ of itself:

$$L_{t1}=D_{t1}+E_{t1}.$$
(8)
The total supply of loans is usually determined by the ratio of asset to equity, so is assumed here as

$$L_{t1}=\frac{1}{\alpha }E_{t1}$$
(9)
with a parameter $`0<\alpha <1`$. The supply of loan to each firm is assumed to be based of collateral of firm asset.

$$L_t^{i(s)}=L_t\frac{K_{t1}^i}{_{j=1}^nK_{t1}^j}$$
(10)
This is where competition between firms is in effect. Then the interest-rate for each firm is determined by $`L_t^{i(s)}=L_t^{i(d)}`$. The bank’s profit is revenue from interests minus the financial cost of the bank

$$\mathrm{\Pi }_t^i=\underset{i=1}{\overset{n}{}}r_t^iL_t^i(\overline{r}_t\omega )L_t$$
(11)
where

$$\overline{r}_t=\underset{i=1}{\overset{n}{}}\frac{L_t^i}{L_{t1}}$$
(12)
and $`\omega `$ is a parameter for profit markup.

The bank faces uncertainty of firms bankruptcy. A firm goes into bankruptcy when its equity becomes negative. Bankrupted firms are replaced with new ones with initial and same balance sheet. The update for the bank’ is thus given by

$$E_t=E_{t1}+\mathrm{\Pi }_t\underset{i=1}{\overset{n}{}}B_{t1}^i$$
(13)
where $`B_{t1}^i`$ is $`A_{t1}^i`$ for every firm $`i`$ bankrupted ($`A_{t1}^i<0`$), otherwise 0.

The dynamics of this system was investigated by simulation in . It was found that the model has, in a stationary phase, skewed distribution of firms sizes with power-law, Laplace distributed growth rates, and that small idiosyncratic shocks generate large aggregate fluctuations. We used the same set of parameters: $`c=1`$, $`\varphi =0.1`$, $`\alpha =0.08`$, $`\omega =0.002`$ except the number of firms $`n=10^5`$ (much larger than $`10^2`$ in original). The stochastic variable $`u_t^i`$ obeys a uniform distribution with support (0,2). The results of simulation are robust with other set of parameters, and with the number of firms.

Our new findings from the simulation is given in Fig. 7. We can observe that the distributions of debt and life-time for bankrupted firms have the same qualitative results as in the real data. These results are also robust for different parameters.

## 4 Summary and Discussion

Firms growth and failure are the two sides of the same coin. In this paper, we have uncovered the phenomenology of firms bankruptcy by using an exhaustive list of Japanese bankruptcies in 1997. Namely, for high-debt regime

1. debt of bankrupted firms follows Zipf distribution
2. debt and size of such firms are related by multiplicative factor independent of the amount of debt
3. bankrupted firms have life-time whose distribution is exponential but correlates with entry rate of new firms

The origin of Zipf law in firms bankruptcy might be well captured in the model , which was shown to reproduce the phenomenological properties in bankruptcy. Indeed the model shares the aspect of self-organized criticality  (see also ) as discussed briefly in what follows.

Failure of a firm can have tremendous influence to economic activities in macroscopic scale. This is because the default of a single firm can propagate its effect through a network of creditor-debtor relationship. In the model, a bankrupted firm leaves the market without paying back its debt and debt commitments to bank. Then the banking system suffers from a capital loss, resulting in the shrink of credit supply. Such shrink accompanied with rise of interest rates is likely to strengthen the competition between firms for seek of more credit and then to deteriorate the firms balance-sheets. But this results in the growing probability of bankruptcy. This is a “domino effect”. The model also shows the opposite case of financially infragile or sound state of firms and bank which leads growth of the total credit and so that of total production (since the system is assumed to have enough money supply when necessary).

During this process of debtor-creditor dynamics, all the firms are competing in their sizes having the risk of bankruptcies. The risk is shared by the bank who attempts to make profit by interest rate. Note that the model has essentially no fine-tuning parameter for having $`\mu =1`$ in eq. (1) and the Zipf law in firms size. In this respect, the model differs from the traditional approach of stochastic processes for firms growth and failure (see  for mid-60’s nice review) including entry-exit processes, collective risk theory, etc. Rather it is the heterogeneous interacting and competing agents that brings about critical state in the aggregate dynamics of firms (see  for this new direction).

We acknowledge Tokyo Shoko Research, Ltd. for providing high-quality data. We are grateful to thank Mauro Gallegati, Gianfranco Giulioni and Nozomi Kichiji for providing information about their model, Hideaki Aoyama and Wataru Souma for discussion, Masanao Aoki for guiding the author to the model, and Katsunori Shimohara for encouragement. The author is indebted to Mauro Gallegati for invaluable comments and suggestions. Supported in part by grants from the Telecommunications Advancement Organization in Japan and from Japan Association for Cultural Exchange.

## Appendix A Bankruptcy data

Bankruptcy or business failure is not a legal term, but should be understood as a commonly used term to describe a state of critical financial insolvency of a debtor. A typical case of bankruptcy is a suspension of banking transactions (see below). What are also termed bankruptcy are: filing to the court for bankruptcy proceedings or a corporate rehabilitation procedure, or informing the aggravation of financial condition of the debtor and delegating creditors to handle it.

Types of bankruptcy are legally divided into legal and private procedures. The former procedure has purpose for rehabilitation or liquidation. Rehabilitation is done by either of Corporate Reorganization Law, Civil Rehabilitation Law (since 2000), Corporate Arrangement under the Commercial Code, or Composition Law (abandoned in 2000). Liquidation is done by Special Liquidation or by “Bankruptcy”. Private procedure, on the other hand, refers to either Suspension of Bank Transactions or Internal Arrangement. The number of each bankruptcy types in 1997 is as follows:

| bankruptcy type | number |
| --- | --- |
| private procedure | |
| Suspension of Bank Transactions | 13850 |
| Internal Arrangement | 501 |
| legal procedure for rehabilitation | |
| Composition Law | 195 |
| Corporate Reorganization | 23 |
| Corporate Arrangement | 12 |
| legal procedure for liquidation | |
| “Bankruptcy” | 1877 |
| Special Liquidation | 68 |
| total | 16526 |

The most frequent case is suspension of bank transactions. Concerning with the bill and check clearings, clearing houses practice a punishment system against the entity who committed dishonored bills or checks in order to maintain the public credibility. The bills or checks which a debtor fails to settle within the term shall be “dishonored”. The second dishonor within six months after the first one shall result in publishing the debtor’s name on the ”the report on suspension of banking transactions.”

Then the debtor shall be barred from transactions with any member financial institutions of the same clearinghouse for two years from the date of suspension, being prohibited from opening trade accounts, drawing bills or checks, and receiving loans. In actuality, this is equivalent to bankruptcy of the firm.

In summary, we regard this case as liquidity problem, and also consider that other legal and private procedures are basically due to the same problem (debt growth).