text:
# Morphology transitions induced by chemotherapy in carcinomas “in situ”

## I Introduction

Cancer is the uncontrolled cellular growth in which neoplastic cells invade adjacent normal tissues and give rise to secondary tumors (metastasis) on tissues or organs distant from its primary origin hanahan . A cancer that remains confined within a messed surface, without break of the underlying basement membrane, is referred to as carcinoma “in situ”. “In situ” carcinoma is characterized by intense cytological atypia, necrosis, and frequent and abnormal mitosis, the tumor cells being arranged in various distinctive patterns mcgee . A malignant tumor is derived from a single cell that, years before the tumor becomes clinically detectable, began an inappropriate pathway of proliferation weinberg . Although cancers are extremely diverse and heterogeneous, a small number of pivotal steps associated to both deregulated cell proliferation and suppressed cell death are required for the development of any and all tumors. Indeed, all neoplasms evolve accordingly to an universal scheme clark ; evan . In the struggle against cancer, surgical removal, chemotherapy and/or radiotherapy are the most commonly used treatment for the complete eradication of the tumor mass. Nowadays, new approaches, such as immunotoxins old , gene navarro , anti-angiogenic folkman  and virus bischoff  therapies, are being developed and have been successfully used in the treatment of several kinds of experimental and human tumors.

Mathematical models are always used as a tentative for describing cancer growth. In particular, numerous models based on classical reaction-diffusion equations have been proposed to investigate the growth of tumor spheroids pettet , cancer progress and its interaction with the immune system bellomo , and the tumor angiogenesis levine ; deangelis . Fractal growth patterns in gliomas (brain tumors) were recently investigated by Sander and Deisboeck sander . Scalerandi et. al., using the local interaction simulation approach (LISA), formulated models for the evolution of avascular tumors scalerandi1  and tumor cords scalerandi2  under nutrient competition. They also analyzed the effects of tumor vascularization scalerandi3  on cancer growth. Recently, we have studied a diffusion-limited model for the growth of carcinomas “in situ” in which cell proliferation, motility and death are locally controlled by growth factors jr1 ; jr2 . This model generates compact, connected and disconnected morphologies characterized by Gompertzian growth in time and distinct scaling laws for the tumor interfaces. In order to generate papillary and ramified morphologies found in many epithelial cancers and trichoblastomas, we were led to analyze the effects of nutrient competition in cancer development jr .

In addition to the vast literature dedicated to tumor growth modeling, many research papers addressing cancer therapies have recently been published. Indeed, cancer cell kinetics under treatments using antimitotic agents kozusco ; montalenti , radiotherapy sachs , virus that replicate selectively in tumor cells wu , anti-angiogenic chemicals scalerandi4 , as wells as the effects of tumor drug resistance and tumor vasculature on chemotherapies jackson , were studied using mathematical models.

In this paper it is analyzed the effects of distinct chemotherapeutic strategies in the model for the growth of avascular tumors proposed by us jr . Specifically, we are interested in possible changes in the tumor patterns, scaling and growth laws reported on the original model triggered by chemotherapies. The paper is organized as follows. In section II the nutrient-limited model for cancer growth and its main results are reviewed. In section III, chemotherapeutic treatments that aim to kill cancer cells are introduced in the framework of the present model and their effects discussed. In section IV, a model for chemotherapy with antimitotic agents is considered. Finally, we draw some conclusions in section V.

## II The nutrient-limited model

Our nutrient-limited model combines macroscopic reaction-diffusion equations, describing the nutrient field concentration, with microscopic stochastic rules governing the actions of individual tumor cells. The basic principles included in the model are cell proliferation, motility and death as well as competition for nutrients among cancer and normal cells. Also, the nutrient concentration field locally controls cell division, migration and death.

### II.1 The model

#### II.1.1 The tissue

The tissue, represented by a square lattice of size $`(L+1)\times (L+1)`$ and lattice constant $`\mathrm{\Delta }`$, is fed by a single capillary vessel at $`x=0`$, the top of the lattice. The capillary is the unique source from which nutrients diffuse through the tissue towards the individual cells. Although a tumor mass is composed of different cell subpopulations clark , the model considers only three types: normal, cancer and tumor dead cells. Any site, with coordinates $`\stackrel{}{x}=(i\mathrm{\Delta },j\mathrm{\Delta })`$, $`i,j=0,1,2,\mathrm{},L`$, is occupied by only one of these cell types. In contrast to the normal cells, one or more cancer cells can pile up in a given site. In turn, dead tumor cells are inert. Thus, each lattice site can be thought of as a group of cells in which the normal, dead and cancer cell populations assume one of the possible values $`\sigma _n(\stackrel{}{x})=\sigma _d(\stackrel{}{x})=0,1`$ and $`\sigma _c(\stackrel{}{x})=0,1,2,\mathrm{}`$, respectively. Accordingly the theory of the monoclonal origin of cancer weinberg , a single cancer cell at $`y=L\mathrm{\Delta }/2`$ and at a distance $`X`$ from the capillary is introduced in the normal tissue. Periodic boundary conditions along the horizontal axis are used. The row $`i=0`$ represents the capillary vessel and the sites with $`i=L+1`$ constitute the external border of the tissue. This geometry is particularly adequate to describe the growth of carcinomas (epithelial tumors) “in situ” because the present model consider that the tumor mass receive nutrients only by diffusion from the capillary vessel. However, the model can be extended to others cancers.

#### II.1.2 The nutrients

The nutrients are divided into two classes: essential to cell proliferation, those necessary to DNA synthesis, and nonessential to cell division. The essential and nonessential nutrients are described by the concentration fields $`N(\stackrel{}{x},t)`$ and $`M(\stackrel{}{x},t)`$, respectively. These nutrient fields obey the dimensionless diffusion equations (see ref. jr  for the complete variable transformations)

$`{\displaystyle \frac{N}{t}}=^2N\alpha ^2N\sigma _n\lambda _N\alpha ^2N\sigma _c`$ (1)
and

$`{\displaystyle \frac{M}{t}}=^2M\alpha ^2M\sigma _n\lambda _M\alpha ^2M\sigma _c`$ (2)
in which differentiated nutrient consumption rates for normal and cancer cells by factors $`\lambda _N`$ and $`\lambda _M`$ are assumed. It is important to notice that the model admits the simplest form for nutrient diffusion, i. e., linear equations with constant coefficients. Also, $`\lambda _N>\lambda _M`$ is used, reflecting the larger affinity of cancer cells for essential nutrients.

The boundary conditions satisfied by the normalized nutrient concentration fields are $`N(x=0)=M(x=0)=1`$, representing the continuous and fixed supply of nutrients provided by the capillary vessel; $`N(y=0)=N(y=L)`$ and $`M(y=0)=M(y=L)`$, corresponding to the periodic boundary conditions along the y-axis; at last, Neumann boundary conditions, $`N(x=L)/x=M(x=L)/x=0`$, are imposed to the border of the tissue.

#### II.1.3 Cell dynamics

In our original model, each tumor cell, randomly selected with equal probability, can carry out one of three actions: division, migration or death. However, in the present work we consider just the accommodation that happens during cell mitosis among the daughter cells. Consequently, each tumor cell can carry out one of two actions.

1. Division. Cancer cells divide by mitosis with probability $`P_{div}`$. If the chosen cell is inside the tumor, its daughter will pile up at that site, and $`\sigma _c(\stackrel{}{x})\sigma _c(\stackrel{}{x})+1`$. Otherwise, if the selected cell is on the tumor border, its daughter cell will occupy at random one of their nearest neighbor sites $`\stackrel{}{x^{}}`$ containing a normal or a necrotic cell and, therefore, $`\sigma _c(\stackrel{}{x^{}})=1`$ and $`\sigma _{n,d}(\stackrel{}{x^{}})=0`$. The mitotic probability $`P_{div}`$ is determined by the concentration per cancer cell of the essential nutrients $`N`$ present on the microenvironment of the selected cell:

$$P_{div}(\stackrel{}{x})=1\mathrm{exp}\left[\left(\frac{N(\stackrel{}{x})}{\sigma _c(\stackrel{}{x})\theta _{div}}\right)^2\right].$$
(3)
The Gaussian term is included in order to produce a sigmoid curve saturated to the unity, and the model parameter $`\theta _{div}`$ controls the shape of this sigmoid.
2. Death. Cancer cells die with probability $`P_{del}`$. Thus, $`\sigma _c(\stackrel{}{x})\sigma _c(\stackrel{}{x})1`$ and $`\sigma _d(\stackrel{}{x})=1`$ when $`\sigma _c`$ vanishes. The cell death probability $`P_{del}`$ is determined by the concentration per cancer cells of the nonessential nutrients $`M`$ present on the microenvironment of the selected cell:

$$P_{del}(\stackrel{}{x})=\mathrm{exp}\left[\left(\frac{M(\stackrel{}{x})}{\sigma _c(\stackrel{}{x})\theta _{del}}\right)^2\right],$$
(4)
i. e., a Gaussian distribution whose variance depends on the model parameter $`\theta _{del}`$.

The biological basis for these cell dynamic rules can be found in reference jr . However, it is worthwhile to notice that from the point of view of the so-called kinetic cellular theory, which provides a general framework for the statistical description of the population dynamics of interacting cells bellomo , the local probabilities $`P_{div}`$ and $`P_{del}`$ can be thought as an effective kinetic cellular model.

#### II.1.4 Computer implementation

The growth model simulations were implemented using the following procedure. At each time step $`T`$, Eqs. (1) and (2) are numerically solved in the stationary state ($`N/t=M/t=0`$) through relaxation methods. Provided the nutrient concentration at any lattice site, $`N_C(T)`$ cancer cells are sequentially selected at random with equal probability. For each one of them, a tentative action (division or death) is randomly chosen with equal probability and the time is incremented by $`\mathrm{\Delta }T=1/N_C(T)`$. The selected cell action will be implemented or not according to the correspondent local probabilities determined by Eqs. (3) or (4). At the end of this sequence of $`N_C(T)`$ tentatives, a new time step begins and the entire procedure (solution of the diffusion equations and application of the cell dynamics) is iterated. The simulations stop if any tumor cell reaches the capillary vessel.

### II.2 Main results

The model reproduces commonly observed tumor morphologies including the papillary, compact and ramified patterns shown in Fig. 2. The nutrient consumption by normal and cancer cells, controlled by the model parameters $`\alpha `$, $`\lambda _N`$ and $`\lambda _M`$, plays a central role in morphology selection. For small values of these parameters, corresponding to growth conditions in which individual cells demand small nutrient supplies, the patterns tend to be compact and circular (Fig. 2a). However, if the mitotic rate of cancer cells is small due to the large amount of nutrients demanded for cell division, generating a significant competition for nutrients, these compact patterns progressively assume papillary-like morphologies (Fig. 2b). At high nutrient consumption rates these papillary patterns become the rule and, for low cancer cell division, continuously transform into thin tips, filaments or chords of cells (Fig. 2c). Also, the model generates patterns with a inner core of died cells for high nutrient consumption or cell division rates (Fig. 2d). As observed in “in vivo” tumors and “in vitro” multicell spheroids pettet , these simulated patterns consist of three distinct regions: a central necrotic core, an inner rim of quiescent cancer cells and a narrow outer shell of proliferating cells.

The tumor patterns generated by our nutrient-limited model were characterized by its gyration radius $`R_g`$, total number of cancer cells $`N_C`$, and number of sites on tumor periphery $`S`$ (including the surface of holes, if any). The gyration radius $`R_g`$ is defined as

$$R_g=\left(\frac{1}{n}\underset{i=1}{\overset{n}{}}r_i^2\right)^{1/2}$$
(5)
where $`n`$ is the number of sites occupied by the pattern (cancer or necrotic cells) and $`r_i`$ is the distance of the occupied site $`i`$ from the tumor mass center. These quantities could be related to clinically important criteria such as progress curves, growth rates (volumetric doubling time) at given radii, proliferative and necrotic fractions of the tumor.

The progress in time of cancer cell populations for all the simulated patterns follows Gompertz curves,

$$N_C(T)=A\mathrm{exp}\left[\mathrm{exp}(k(TT_c))\right],$$
(6)
representing an initial exponential growth that is saturated in larger times.

Alternatively, as a function of the number of sites occupied by the pattern $`n`$, both $`R_g`$ and $`S`$ obey power law scalings given by $`R_gn^\nu `$ and $`Sn^\sigma `$, respectively. For solid patterns these exponents are $`\nu 0.5`$ and $`\sigma 0.5`$, corresponding to effective circular and non fractal patterns. As the nutrient consumption increases, the patterns tend to papillary-like shapes for which the exponent $`\sigma `$ increases towards the value $`1`$ and the exponent $`\nu `$ varies in the range $`[0.50,0.60]`$, indicating a fractal morphology for the tumor.

## III Cytotoxic chemotherapeutic drugs

The primary aim of the antitumoral treatment with chemotherapeutic drugs is to kill or at least to stop the proliferation of the cancer cells. In general, the drugs should act only on proliferating cells, mainly the cancer ones. However, drugs also destroy proliferating normal cells promoting several collateral effects hellman . Indeed, epithelial cells from the respiratory and gastro-intestinal systems, which frequently reproduce in order to substitute their dead counterparts, are strongly affected. In this section, we analyze a simple chemotherapeutic model in which the complex details of the cell-cycling responses to the drugs are taking into account by an effective kinetic model.

### III.1 The model

As used in previous works panetta2 , the chemotherapy is modelled by a periodic delivery of cytotoxic drugs through the same capillary vessel supplying the nutrients to the tissue. Several cytotoxic drugs and their properties were exhaustively studied such as amsacrine, cisplatin, cyclophosphamide, cytarabine, mustine, and anthracycline oxford . Here, the numerous barriers involved in tumor drug delivery jain  were not considered, and the treatment begins when the tumor mass contains $`N_0`$ cancer cells. When a dose is applied, the drug concentration in the capillary assumes a value $`\stackrel{~}{Q}_0`$. By hypothesis, this concentration level decays exponentially in time, simulating the gradual drug elimination by the organism. New doses are periodically administered at time intervals $`\tau `$. So, the drug concentration at the capillary, $`\stackrel{~}{Q}(T)`$, is a function of time written as

$$\stackrel{~}{Q}(T)=\{\begin{array}{c}0,\text{ if }T<T_0\hfill \\ \stackrel{~}{Q}_0\mathrm{exp}[(Tl\tau )/T_\times ],\text{ if }T_0+l\tau T<\hfill \\ T_0+(l+1)\tau \hfill \end{array},$$
(7)
where $`l=0,1,2\mathrm{}`$ is the number of applied doses, $`T_\times `$ is the characteristic time for the drug elimination by the organism and $`T_0`$ is the time at which the treatment begins. The functional form of $`\stackrel{~}{Q}(T)`$ is shown in Fig. 1.

As the nutrients, drugs diffuse from the capillary vessel toward the individual cells and their concentration $`Q(\stackrel{}{x},t)`$ at any lattice site at each time step is given by the stationary solution of the diffusion equation

$$\frac{Q}{t}=D_Q^2Q\mathrm{\Gamma }^2Q\lambda _Q\mathrm{\Gamma }^2Q\sigma _c.$$
(8)
The first, second and third terms on the right hand side represent drug diffusion, natural degradation and absorption by cancer cells, respectively. Again, this equation is the simplest one describing the diffusion phenomenon. Eq. (8) is written using the same dimensionless variables of (1) and (2). Also, the parameters $`\stackrel{~}{Q}_0`$ and $`D_Q`$ can be make equal to the unity without generality loss. The boundary conditions are the same used for the nutrient fields, except at the capillary vessel where the concentration at each time step is given by Eq. (7).

Finally, a single change is introduced into the cell dynamics described in section II: an additional chance of death occurs whenever a cancer cell enters mitosis. Thus, every time a cancer cell divides, each one of the generated cells can die with probability

$$P_{del}^{(Q)}(\stackrel{}{x})=1\mathrm{exp}\left[\left(\frac{Q(\stackrel{}{x})}{\sigma _c(\stackrel{}{x})\theta _{del}^{(Q)}}\right)^2\right].$$
(9)
The parameter $`\theta _{del}^{(Q)}`$ controls the cell sensitivity to the drugs. Also, since in the original nutrient-limited model the normal cells do not take part on the cell dynamics, i. e., they do not divide or die, we disregard the chemotherapeutic effects on normal cells.

### III.2 Results

The main aim of this work was to investigate the effects of treatments on the various morphologies, scaling and growth laws observed in the original model. Thus, for each morphology, determined by the fixed parameters reported in table 1, the treatment parameters, namely, $`\tau `$ (the dose period), $`\theta _{del}^{(Q)}`$ (drug efficiency) and $`N_0`$ (tumor size at the treatment beginning) were varied. Such parameters can be directly tested in the laboratory. The remaining parameters ($`\lambda _Q`$, $`\mathrm{\Gamma }`$ and $`T_\times `$) associated to drug diffusion and elimination were also fixed in all the simulations.

In Fig. 2, the corresponding patterns for compact, papillary, and ramified morphologies with and without treatment are shown. In this figure the treatment is not able to cease tumor development. As can be seen, the morphological tumor patterns do not change under mild treatment. However, the regions occupied by the tumors are larger than those without treatment. These results suggest that the direct attack to the tumor might be an inadequate treatment strategy. Indeed, the more invasive is the tumor, higher is the possibility of the cancer cells to reach the capillary vessel and, consequently, metastasize successfully. Such result is in agreement with the Israel’s claim lucien  that cancer cells trigger adaptation mechanisms in stress circumstances similar to those observed in bacterial colonies. Moreover, in our model neither genetic or epigenetic changes are necessary to explain the increase of tumor aggressiveness. It naturally emerges from the growth rules. A tumor submitted to successive chemotherapeutic treatments which do not lead to its complete eradication may progressively become more resistant, aggressive and malignant.

Very similar scaling laws $`R_gn^\nu `$ and $`Sn^\sigma `$ for the treated tumors were observed ($`\nu `$ and $`\sigma `$ values are reported in jr ). The small differences in the exponents values for the number of tumor peripheral sites vanishes at the asymptotical limit of tumor size. This exponents invariance suggests that the fractal morphology is a robust aspect of these tumors and cannot be changed by small perturbations in the cell microenvironment. Obviously, the scaling laws for $`R_g`$ and $`S`$ are not defined for the tumors that cease their growth.

We have also studied the influence of the treatment parameters ($`\tau `$, $`\theta _{del}^{(Q)}`$, $`N_0`$) on the tumor growth curves. In Fig. 3 the evolution in time of the total number of cancer cells for compact tumors is shown. Depending on dose period $`\tau `$ the tumors may disappear, saturate their growth or progress according to power laws. Actually, the time interval between two consecutive doses is a fundamental clinical feature determining the treatment success. The power law regime observed when $`\tau `$ is large means a slow progress, neatly contrasting with the exponential growth present in the Gompertz law describing tumor progress without chemotherapy. The tumor size at the beginning of the treatment is another important factor. Indeed, the exponents of the power laws characterizing the tumor growth are smaller when the initial tumor size is bigger. This scenario seems to be universal. Sometimes, as shown Fig. 3, compact tumors which began to receive drug doses at regular intervals of $`\tau =5`$ when it contained $`N_0=10^4`$ cancer cells were eliminated slower than another one that began to be treated with $`N_0=5\times 10^4`$ cancer cells. However, this is not the rule. For example, simulations of papillary tumors indicate that smaller tumors are faster eliminated when $`\tau =5`$. Finally, since the growth law exponents depend on the parameter sets used, they are not universal. Concerning the $`\theta _{del}^{(Q)}`$ parameter, it just modifies the treatment efficiency.

Fractal tumors are more resistant to treatments. In Fig. 4, the growth curves $`N_C\times T`$ are drawn for compact, papillary, ramified and necrotic morphologies shown in Fig. 2. One can see that the more fractal is the tumor larger is the required time to eliminate it. Indeed, the lower growth rates of fractal tumors imply in a large fraction of cancer cells maintained at a quiescent state. So, since chemotherapeutic drugs considered in this model only act in the dividing cells, the treatment becomes inefficient when the major fraction of the cancer cells are quiescent.

In addition, we have studied tumor patterns exhibiting a central necrotic core. The results for cancer growth and the correspondent power laws are similar to those found for compact patterns. The density of cancer cells and their average division rates through the growth patterns are not significatively altered when the tumors are submitted to mild treatments (long time intervals between consecutive drug doses). Moreover, their growth patterns are very similar to those found in the untreated counterparts. Finally, tumors exhibiting necrotic cores are more easily eliminated when shorter periods of the drug administration are considered (Fig. 4).

## IV Antimitotic chemotherapeutic drugs

Several drugs used in cancer chemotherapies do not kill cancer cells. Instead, they aim to stop the cell cycle in specific checkpoints. As a consequence of that, the tumors cease or slow down their growth rate. Examples of such drugs include the antimitotic agent Curacin A that blocks the cell cycle in the M-phase, mitomycin C, doxorubicin, and aclarubicin, among others oxford . In order to analyze the effects of drugs that inhibit cell division on tumor patterns, we introduce a very simple change in the cell dynamics of the model described in the previous section. Instead of an additional death probability given by Eq. (9), the division probability is modified by including an additional term in Eq. (3). So, the new division probability becomes

$$P_{div}=1\mathrm{exp}\left[\left(\frac{N}{\sigma _c\theta _{div}}\right)^2+\left(\frac{Q}{\sigma _c\theta _{div}^{(Q)}}\right)^2\right],$$
(10)
in which the first and second terms in the exponential argument compete between themselves stimulating and inhibiting the cell mitosis, respectively. The parameter $`\theta _{div}^{(Q)}`$ controls the drug influence on the cell division. Obviously, if the exponential argument is larger than the unity then $`P_{div}0`$.

The biological interpretation of the present model is very distinct from that considered in the last section. Now, the drugs modify the intrinsic characteristics of the cancer cells and, consequently, the tumors must behave differently from those treated with cytotoxic agents. Indeed, as one can see in Fig. 5, a morphological transition for the tumor patterns when the dose periods decrease is observed . If the period is sufficiently short the tumor size remains constant during the therapy. Obviously, the critical period depend on other model parameters, specially on $`\theta _{div}^{(Q)}`$. In order to see this dependence, was evaluated the critical value of $`\theta _{div}^{(Q)}`$ for which the tumors cease their growth with probability $`1`$ as a function of $`\tau `$. A mean field analysis of Eq. (10) provides $`\theta _{div}^{(Q)}\mathrm{exp}(\tau /T_\times )`$, in agreement with the simulations as indicated in Fig. 6. Thus, it was found an exponential decay with an universal characteristic length $`T_\times `$ that is independent of the other model parameters. This law, relating two important clinical parameters, is valid for an wide set of parameter values. In addition to the compact tumors shown in figure 5, it was also simulated papillary tumors under antimitotic treatment, and a similar morphology transition was observed.

Fig. 7 shows that these morphological transitions occurred at well defined $`\theta _{div}^{(Q)}`$ values. Near the critical value, instabilities in the reaction-diffusion equations lead to a branching growth sander . Below this threshold the cell death rate equals cell division, and the tumor growth effectively stops. This qualitative behavior near the transition will probably not be observed in experimental tumors, since a very reliable determination of $`\theta _{div}^{(Q)}`$ is very difficult. On the other hand, experiments exhibiting both regimes, above and below the critical value, are realizable today, due to the great advances in drug delivery control. Such assays might corroborate the complex behaviors predicted by our model.

The antimitotic treatment affects both the growth and scaling laws. The compact patterns in the original model, for which the gyration radius and number of occupied sites on the tumor border scale as the square root of the total number of occupied sites (i. e. $`\nu =\sigma =1/2`$), become progressively more fractal as the dose intervals are decreased. Consequently, it was obtained $`\nu >1/2`$ and $`\sigma >1/2`$ indicating fractal tumors with a fractal dimension given by $`d_f=1/\nu `$. In general, all the tumors become more fractal when submitted to antimitotic treatments. So, there is a neat contrast with the invariance of the fractal dimension of the growing tumors under the cytotoxic treatment described in the previous section. In Fig. 8 the increase in the number of cancer cells relative to the treatment beginning ($`N_CN_0`$) is plotted, for distinct dose intervals, as a function of the time afterwards the treatment beginning ($`TT_0`$). These curves suggest power laws for the time evolution of the number of cancer cells with a weak dependence on the $`\tau `$ parameter. Indeed, we have $`(N_CN_0)(TT_0)^\beta `$ whit $`\beta [1.2,1.4]`$. As $`\tau `$ decreases, the number of tumor cells initially decays for a while, but subsequently recovers its growth (inset of Fig. 8). Below a given dose interval, $`N_C`$ decreases monotonically up to a constant value.

In order to analyze the cell division rates through the tumor, we computed the average cancer cell density ($`\sigma _c`$) and cell division rate ($`w_{div}`$) along a longitudinal cut across the growth pattern. These average values were plotted as a function of the distance from the capillary vessel. In Fig. 9 these plots are shown for compact patterns treated with distinct $`\tau `$ values. For larger $`\tau `$ values (Fig. 9(a)) both, division rate and cancer cell density, have sharp maxima on the tumor border in front and opposite to the capillary vessel. Notice that the peaks for division rate are significatively narrower than those for cancer cell density, demonstrating that the proliferative fraction of the tumor comprises just a small rim located at the tumor border. As the treatment is intensified (shorter $`\tau `$ values are considered) the cell densities through the tumors become more uniform, whereas the division rates continue exhibiting sharp maxima at the tumor borders (Figs. 9(b) and (c)). If $`\tau `$ is sufficiently short in order to halt the tumor growth, the division rates and cancer densities become uniform along the tumor (Fig. 9(d)). So, in these cases the division rates are counterbalanced by the death rates, leading to a vanishing net rate of cancer growth. All these results show that the tumor patterns are altered when the mitotic properties of the cancer cells are modified by external agents.

## V Conclusions

In the present work, a reaction-diffusion model to simulate the effects of chemotherapies on the growth of carcinomas “in situ” have been studied. The model includes cell death and division, competition among cancer and normal cells by nutrients and periodical drug administration. Two kinds of chemotherapies, using cytotoxic and antimitotic drugs, were modelled.

In the cytotoxic model the tumors can be completely eradicated, cease their growth or grow continuously. The last case occurs when the treatment is inefficient or the intervals between consecutive doses of drugs are large. Moreover, the morphologies and scaling laws of the growing tumor patterns are preserved. In contrast, for therapeutic approaches using antimitotic agents a morphological transition in the tumor patterns was observed. The growth patterns become progressively more fractal as more effective treatments (shorter intervals between consecutive doses and/or more efficient drugs) are considered.

Such morphological transitions are similar to recent studies claiming that bacterial colonies exposed to non-lethal concentrations of antibiotics exhibit drastic changes in their growth patterns ben-jacob . For bacteria, these changes were imputed to variations in bacterial properties such as metabolic load and chemotaxis. In turn, normal and cancer cells cultured in monolayer and collagen gel exhibit a dynamical transition in their aggregation regimes as an adaptive response to the growth constraints imposed by high cell population density or long permanence in culture mendes . Again, the results obtained in the present work are in agreement with the point of view that cancer cells can develop an integrated defense program against stress situations similar to the response of bacterial colonies facing severe and sustained threats lucien . However, as far as we are concerned, there are no reports on cancer literature concerning transitions morphologies in the histological patterns of tumors submitted to chemotherapy. Thus, formal models like the one proposed in this paper, generally not familiar for the most of the biomedical researchers complexity , can guide and refine new experiments intended to analyze such morphological transitions. Currently, experiments “in vitro” with cancer cells are being performed in our laboratories in order to investigated those morphological transitions.

It is important to mention that, usually, real effective treatments employ several therapeutic methods simultaneously hellman . So, a more realistic chemotherapeutic model should consider combined cytotoxic and antimitotic treatments. Nevertheless, in these cases relevant features such as the robustness of the tumors fractal scaling and the morphology transitions can be masked. At last, we are modelling others cancer therapeutic strategies, such as virus and immunotoxins therapies, as well as the combination of distinct treatments by using reaction-diffusion models like that considered in this paper.

###### Acknowledgements.

S.C. Ferreira Jr. would like to thank the UFV Physics Department (Viçosa, Brazil) for the welcoming hospitality and the Professor J. G. Moreira (Departamento de Física, UFMG, Brazil) for fruitful discussions. This work was partially supported by the CNPq an FAPEMIG -Brazilian agencies.


format: text

text:
# Alkali-metal ion in rare gas clusters: global minima

## 1 Introduction

Rare gas clusters are important and often studied aggregates of atoms . Their mass spectra show the presence of magic numbers which correspond to the most stable geometrical structures . Because of the simple form of van der Waals chemical bonding rare gas clusters are very suitable for theoretical studies. For each cluster size different isomeric structures exist. They are local minima in the multidimensional potential energy surface. The number of these isomers increases rapidly with the size of a cluster. Several methods of a search for global minima were invented . When the size of clusters increases it becomes more difficult to find the lowest minimum on the potential energy surface. Optimization methods locate only the candidate for a global minimum.

The problem of interaction between various clusters and impurities generated intensive research . For example, impurities in helium nanoparticles were studied . One goal was to determine whether impurity ions, atoms, or molecules preferably reside in the interior of the cluster, or on its surface. In general, the impurities substantially perturb the cluster matrix. It was found that a single metal atom drastically modifies the properties of silicon clusters . Therefore, it is interesting to explore the shape deformation and the energy change induced by impurities. In this context rare gas clusters are especially interesting. They are inert and could play a role of the nanoscale chemical laboratory. A global optimization is one of theoretical methods suitable for investigation of the impurity-cluster interaction.

The most successful methods of global optimization for clusters are genetic algorithms  and a basin-hopping method . The genetic algorithm methods select a fraction of the initial population of relaxed cluster structures as “parents” using their energies as the criteria of fitness. The next generation of cluster structures (“children”) is formed by “mating” these parents and the process is repeated until the minima are found. The basin-hopping approach of global optimization for clusters was developed by Wales and Doye . This method uses the energy minimization and the Monte Carlo simulation. The potential energy is transformed to a new surface which has the form of a multi-dimensional staircase. The steps correspond to the basins of attraction surrounding energy minima . Each basin of attraction contains the set of configurations which after geometry optimization generate a particular minimum. The basin-hopping method is very powerful and, for example, finds all known energy minima for the Lennard–Jones clusters containing up to $`110`$ atoms . A similar method was used in the study of protein folding . The Lennard–Jones nanoparticles play a role of the test system for global optimization methods. Molecule-doped rare gas clusters, in particular $`H_2`$ in $`Ar_N`$ , $`Cl_2`$ in $`Ar_N`$ , and $`NO`$ in $`Ar_N`$ , were also studied. In general, the clusters with impurities are less studied by optimization methods than homonuclear ones. A small perturbation, for example one impurity atom, or ion, in a cluster may shed light on the global minima formation. In addition, new morphologies in the set of cluster structures may appear for some impurities. The methods of global optimization produce the candidates for minima of the free energy only close to $`T=0`$ K. At higher temperatures kinetics effects are important and may change the order of stability of cluster isomers.

In this work I present the results of structural optimization for model heteronuclear clusters consisting of one impurity ion and of up to $`N=79`$ rare gas atoms. The Lennard–Jones and Mason–Schamp potentials are used to describe interactions. The basin-hopping method is applied for structural optimization. The results show that an impurity ion is at the surface for a small number of rare gas atoms. It is always trapped in a cluster for $`N8`$ atoms. As in homonuclear Lennard–Jones clusters, several particularly stable structures are found. In the following the computational method is described in Sec. II. Results and discussion are presented in Sec. III. A summary and conclusions are given in Sec. IV.

## 2 Computational method

The pairwise additivity and the spherical symmetry of the interaction for rare gas atoms are reasonably well described by the Lennard–Jones potential. This potential is one of the most used models of interaction. It is given by

$$V(r)=4ϵ\left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right],$$
(1)
where $`ϵ`$ is the depth of the potential energy minimum, and $`2^{1/6}\sigma `$ is the equilibrium pair separation . The $``$ $`r^6`$ is the leading term in the dispersion energy interaction for neutral atoms. Mason and Schamp proposed the additional term $``$$`r^4`$ to model the interaction of an ion and a neutral rare gas atom . The Mason–Schamp potential is given by

$$V(r)=\frac{ϵ}{2}\left[(1+\gamma )\left(\frac{\sigma }{r}\right)^{12}4\gamma \left(\frac{\sigma }{r}\right)^63(1\gamma )\left(\frac{\sigma }{r}\right)^4\right].$$
(2)
The coefficients in this potential were determined from the ion mobility measurements in rare gases and calculations within the kinetic theory. Three parameters in algebraic form (2) specify the depth and position of the potential energy minimum, and the relative strength of various terms. The Mason–Schamp potential was used in the studies of impurities in rare gas crystals . For example, the parameters of this potential for the $`Cs^+`$$`Xe`$ system are: $`ϵ=0.106`$ eV, $`\sigma =3.88\AA `$, and $`\gamma =0.2`$. Figure 1 shows the Mason–Schamp and Lennard–Jones potentials. The Mason–Schamp potential was taken here in its reduced form as a model of a class of the alkali-metal ion and rare gas atom interaction. In the same sense the Lennard–Jones potential represents interaction between rare gas atoms in investigations of the global minima for homonuclear clusters . The Lennard-Jones and Mason–Schamp potentials in the inset of Fig. 1 are drawn in reduced units where $`ϵ`$ and $`\sigma `$ are chosen as the units of energy and distance. In the calculations the total potential was modelled by a pairwise sum of potentials given by Eq. (1) for each rare gas/rare gas pair of atoms, and by Eq. (2) for each rare gas/ion pair of particles. For both potentials reduced units for $`ϵ`$ and $`\sigma `$ were employed.

I used a basin-hopping method and a modified program gmin of Wales and collaborators . The constant temperature Monte Carlo method based on the standard Metropolis scheme  is applied. Energy minimization is performed using the limited memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm for the unconstrained minimization problem . For each cluster size five series of runs were done starting from different configurations. Each run consisted of $`30000`$ Monte Carlo steps at a constant reduced temperature of $`0.8`$. Control runs for homonuclear Lennard–Jones clusters have been carried out, and the same global minima were found as by Wales and Doye . For sizes $`N<75`$, direct optimizations from random configurations were able to produce known minima for homonuclear Lennard–Jones clusters within less than $`10000`$ Monte Carlo steps. For $`N75`$, additional effort was necessary, either the increase of the number of steps, and/or seeding with configurations derived from the best structures around a chosen cluster size. For both homonuclear and heteronuclear clusters the lowest-energy structures were found at different Monte Carlo steps in five series of runs. Sometimes, different lowest-energy structures were found in these runs. Differences between various runs were especially pronounced for $`N75`$. Similar type of potentials are used for the atom-atom and atom-ion interactions and it is plausible that the same procedure is adequate for unknown energy minima of heteronuclear clusters. However, in heteronuclear systems the forces are determined not only by distances between particles, but their chemical identity is also important. In homonuclear clusters swapping the position of the particles does not change structural properties of the aggregate. For the system and the method presented in this work (where only one impurity particle, two similar potentials, and the Monte Carlo minimization method are used) the optimization tasks for homonuclear and heteronuclear clusters are analogous. In theoretical analysis of global optimization and NP-hard problems, the case of heteronuclear clusters is more difficult .

## 3 Results and discussion

Figures 2 and 3 show several configurations of an alkali-metal impurity in a rare gas cluster.

This visualization is performed using the Rasmol package . It is found that all lowest-energy structures of heteronuclear clusters are compact. The spacial distributions of particles does not change much in comparison with homonuclear Lennard–Jones clusters. The ion-atom distance is smaller than the atom-atom one. The lengths of the bonds and the angles between them change slightly. The ratio of the numbers of the atom-atom and ion-atom bonds increases with the cluster size. Therefore, a difference between corresponding homonuclear and heteronuclear clusters decreases when their size increases. The energy minima are presented in Table 1.

The situation where the alkali-metal ion is at the surface of a cluster is shown in Fig. 2(a-f). Figures 2(f) and 2(g) present configurations around the transitional regime of the ion. In Fig. 2(f) the impurity ion is deeply buried in the surface of the cluster, and lies closer to its centre than one rare gas atom. The ion is in the interior of the cluster the first time for eight rare gas atoms, as shown in Fig. 2(g). The impurity ion is also always in the cluster for more than eight rare gas atoms (Figs. 2(h-j) and Fig. 3). When the alkali-metal ion is solvated by rare gas atoms it is not always situated in the centre of the cluster. For example, the impurity ion takes a non-central position in clusters shown in Figs. 3(b-d).

The energies per particle are presented in Fig. 4(a). This quantity varies rather monotonically. However, mild minima occur for certain sizes and show particularly stable structures. This is further confirmed by the first energy differences $`\mathrm{\Delta }_1=E(N)E(N1)`$ shown in Fig. 4(b). Two minima at $`N=13`$ and $`N=55`$ correspond to well known filled icosahedral shells. In addition, minima occur for $`N=19,39,46,49,70,79`$. These most stable structures of rare gas clusters with an alkali-metal ion are shown in Fig. 3. Minima at $`N=19,46,49,70,79`$ also exist for Lennard–Jones clusters (for example, see Fig. 4 in ). However, in homonuclear Lennard–Jones clusters minimum exist for $`N=38`$. Figure 4(b) shows that for heteronuclear clusters corresponding minimum is $`N=39`$. Shapes of $`N=38`$ and $`N=39`$ (see Fig. 3(c)) heteronuclear clusters are similar to each other and to Lennard–Jones $`N=39`$. The Lennard–Jones $`N=38`$ cluster is slightly more compact (see Fig. 1 in ). A shape of the $`N=39`$ structure shown in Fig. 3(c) is closer to the third lowest energy minimum of Lennard–Jones $`N=38`$, than to the first and second one. In connection with the complex double-funnel energy landscape of the non-icosahedral $`38`$-atom Lennard–Jones cluster , it is interesting to point out that this magic size in heteronuclear nanoparticles is replaced by $`N=39`$. Other non-icosahedral Lennard–Jones global minima occur for $`N=75`$$`77`$, and these clusters are Marks decahedra . For clusters with the impurity ion, sizes $`N=75`$$`77`$ are decahedral in origin. As for icosahedral minima, these structures are distorted in comparison with homonuclear Lennard–Jones clusters. The flat region in Fig. 4(b) exist for $`51N55`$. In these clusters the impurity ion is situated in the central position and the magic icosahedral $`N=55`$ structure builds up. For homonuclear Lennard–Jones clusters a such flat region is $`52N55`$ , and it also corresponds to the filling of the magic $`N=55`$ structure.

## 4 Conclusions

The global optimization based on a basin-hopping method is carried out to study structural properties of model nanoparticles consisting of one impurity alkali ion and a rare gas cluster. The catalog of energy minimum candidates for these clusters containing up to $`80`$ particles is presented. It is found that the impurity ion is lying in the surface of a nanoparticle for a small number of atoms. An alkali-metal ion is solvated by rare gas atoms when their number is greater or equal than eight. The predicted magic numbers $`13,19,46,49,55,70,79`$ of heteronuclear clusters are the same as for homonuclear Lennard–Jones ones. The magic structure $`N=39`$ of heteronuclear clusters replaces $`N=38`$ of homonuclear nanoparticles. Therefore a small perturbation, such as the presence of the closest alkali-metal ion neighbour from the periodic table, may change the list of magic sizes for rare gas clusters. The system studied in this work is useful as a test example for theoretical analysis of global optimization and NP-hard problems for heteronuclear clusters . The methods of global optimization are helpful as a guide in low-temperature experimental studies of nanoparticles with impurities.

This work has been carried under the HR-MZT project 0119255 “Dynamical Properties and Spectroscopy of Surfaces and Nanostructures”.


text:
# Hadrons with charmed quarks in matter

## 1 Introduction.

In dense and hot nuclear matter the light quark condensates $`q\overline{q}`$ might be substantially reduced. This affects the light quark content of mesons and baryons and therefore results in medium modification of hadron properties Hatsuda ; Weise ; Tsushima1 . Even if the changes in quark condensates are small, the absolute difference between the in-medium and bare masses of hadrons is expected Tsushima1 ; Sibirtsev1  to be larger for heavier hadrons.

The charmed mesons, which consist of light $`q,\overline{q}`$ and heavy $`c,\overline{c}`$ quarks, are considered suitable probes of in-medium modification of hadron properties. Similaly to $`\overline{K}`$ ($`\overline{q}s`$) and $`K`$ ($`q\overline{s}`$) mesons, the $`D`$ ($`\overline{q}c`$) and $`\overline{D}`$ ($`q\overline{c}`$) satisfy different dispersion relations in matter because of the different sign of $`q`$ and $`\overline{q}`$ vector couplings Tsushima1 . While the $`D`$ meson mass is reduced in nuclear matter, the $`\overline{D}`$ mass is raised, as is illustrated by Fig.1a). Calculations with the Quark Meson Coupling model (QMC) Guichon  show that already at normal nuclear density $`\rho _0`$ the mass splitting between $`D`$ and $`\overline{D}`$ mesons is about 160 MeV.

It was proposed Sibirtsev1  that the modification of the $`D`$ meson in nuclear matter can be identified by enhanced subthreshold production of open charm in $`\overline{p}A`$ annihilation. Because of charm conservation, $`D`$ and $`\overline{D}`$ mesons are produced pairwise. As is shown in Fig.1b) the sum of $`D`$ and $`\overline{D}`$ masses depends substantially on nuclear density. The downward shift of the $`D\overline{D}`$ threshold at $`\rho _0`$ is $``$–100 MeV. A QCD sum rule estimate Hayashigaki  predicts about the same shift.

Furthermore, an attractive $`D`$-nucleus potential can be measured by investigating charmed mesic nuclei Tsushima1 . The reduction of $`D`$ mass in matter might affect open charm production Cassing  and $`J/\mathrm{\Psi }`$ suppression Sibirtsev2  in relativistic heavy ion collisions.

In contrast to open charm, charmonium mesons consist of heavy $`c\overline{c}`$ quarks and can be affected only by gluonic condensates. It was expected that the modification of properties of heavy quarkonia in matter would be almost negligible Brodsky . In that case, the overall $`D\overline{D}`$ mass at some nuclear density might cross the masses of excited charmonium states, as is illustrated by Fig.1b).

The crossing of the charmonium states levels and the $`D\overline{D}`$ threshold in nuclear medium might result in the melting Digal  of excited charmonium mesons.

Here we investigate the modification of widths of charmonium states at finite nuclear density with $`D`$ and $`\overline{D}`$ in-medium masses predicted by QMC Tsushima1 . We consider simultaneous modification of both $`D`$ and $`\overline{D}`$ mesons, as is illustrated in Fig.1 and study $`\mathrm{\Psi }^{\prime \prime }`$(3770), $`\mathrm{\Psi }^{}`$(3686) and $`\chi _{c2}`$ decay into $`D\overline{D}`$ in nuclear matter. The $`\chi _{c0}`$ modification is not discussed since the $`D\overline{D}`$ threshold does not cross its mass even at $`\rho `$=3$`\rho _0`$, as is shown in Fig.1. The $`\chi _{c1}D\overline{D}`$ decay is suppressed by parity conservation.

## 2 $`\mathrm{\Psi }^{\prime \prime }`$(3770).

The $`\mathrm{\Psi }^{\prime \prime }`$ charmonium lies above the $`D\overline{D}`$ threshold in free space and its dominant decay width into $`D\overline{D}`$ channel is given by

$$\mathrm{\Gamma }_{\mathrm{\Psi }^{\prime \prime }D\overline{D}}=\frac{g_{\mathrm{\Psi }^{\prime \prime }D\overline{D}}^2}{3\pi }\frac{q^3}{m_\mathrm{\Psi }^2},$$
(1)
where $`m_\mathrm{\Psi }`$ is the $`\mathrm{\Psi }^{\prime \prime }`$(3770) mass and $`q`$ is the $`D`$ meson momentum in the charmonium rest frame,

$$q=\frac{[(m_\psi ^2m_D^2m_{\overline{D}}^2)^24m_D^2m_{\overline{D}}^2]^{1/2}}{2m_\mathrm{\Psi }},$$
(2)
with $`m_D`$ and $`m_{\overline{D}}`$ the masses of $`D`$ and $`\overline{D}`$ mesons, respectively, while the coupling constant $`g_{\mathrm{\Psi }^{\prime \prime }D\overline{D}}`$=14.89 is fixed by the vacuum decay width $`\mathrm{\Gamma }_{\mathrm{\Psi }^{\prime \prime }D\overline{D}}`$=23.6 MeV.

If $`\mathrm{\Psi }^{\prime \prime }D\overline{D}`$ coupling does not change in matter the modification of the $`\mathrm{\Psi }`$(3770) width is entirely given by $`D`$ and $`\overline{D}`$ in-medium masses and is determined by the phase space dependence of the charmonium decay width. In that case the in-medium $`\mathrm{\Psi }`$(3770) width depends substantially on nuclear matter density, as is shown by the solid line in Fig.2a).

Within the $`{}_{}{}^{3}P_{0}^{}`$ model Busetto  the $`\mathrm{\Psi }^{\prime \prime }D\overline{D}`$ coupling itself depends on the $`D`$ and $`\overline{D}`$ masses via

$`g_{\mathrm{\Psi }^{\prime \prime }D\overline{D}}^2={\displaystyle \frac{\pi ^{3/2}\mathrm{\hspace{0.17em}2}^{11}\gamma ^2m_\mathrm{\Psi }}{5\beta ^7\mathrm{\hspace{0.17em}3}^{10}}}[(q^2+m_D^2)(q^2+m_{\overline{D}}^2)]^{1/2}`$
$`\times (15\beta ^22q^2)^2\mathrm{exp}(q^2/6\beta ^2),`$ (3)
where the oscillator length scale $`\beta `$=360 MeV is fixed by light mesons decays Busetto , while the interaction strength $`\gamma `$=0.33 is determined by $`\mathrm{\Psi }^{\prime \prime }D\overline{D}`$ decay.

Finally, the dashed line in Fig.2 shows the dependence of the $`\mathrm{\Psi }^{\prime \prime }D\overline{D}`$ decay width on nuclear matter density resulting from the $`{}_{}{}^{3}P_{0}^{}`$ model. At normal nuclear density $`\rho _0`$= 0.16 fm-3 the $`\mathrm{\Psi }^{\prime \prime }(3770)`$ in-medium width almost saturates $`\mathrm{\Gamma }_{\mathrm{\Psi }^{\prime \prime }D\overline{D}}`$90 MeV, which is $``$3.8 times larger than in vacuum. Furthermore, the $`{}_{}{}^{3}P_{0}^{}`$ result substantially differs from the phase space estimate.

The $`\mathrm{\Psi }^{\prime \prime }`$ modification in nuclear matter might be studied by dilepton spectroscopy from $`AA`$ as well as $`\overline{p}A`$ interactions, since the effect is measurable already at normal nuclear densities. The $`\mathrm{\Psi }^{\prime \prime }(3770)`$ charmonium does not decay into $`J/\mathrm{\Psi }`$ and thus its modification in matter can not be considered as an additional source of $`J/\mathrm{\Psi }`$ suppression in heavy ion collisions.

## 3 $`\mathrm{\Psi }^{}`$(3686).

The $`\mathrm{\Psi }^{}`$ charmonium lies below the $`D\overline{D}`$ threshold in free space, but the $`D\overline{D}`$ decay channel becomes open at nuclear density $`\rho `$0.05 fm-3, as illustrated by Fig.1. The $`\mathrm{\Psi }^{}`$ is narrow, its total width is 0.277 MeV and partial decay into $`J/\mathrm{\Psi }`$ accounts for $``$54% of the total width.

Although, the $`\mathrm{\Psi }`$(3686) coupling to $`D\overline{D}`$ is not directly accessible, it can be estimated within the framework of Vector Meson Dominance model as

$$g_{\mathrm{\Psi }^{}D\overline{D}}^2=\frac{16\pi \alpha ^2}{27}\frac{m_\mathrm{\Psi }}{\mathrm{\Gamma }_{\mathrm{\Psi }^{}e^+e^{}}},$$
(4)
where $`\alpha `$ is the fine structure constant, $`m_\mathrm{\Psi }`$ is the $`\mathrm{\Psi }^{}`$(3686) mass and $`\mathrm{\Gamma }_{\mathrm{\Psi }^{}e^+e^{}}`$=2.35 KeV is the radiative $`\mathrm{\Psi }^{}e^+e^{}`$ decay width. Finally, the $`\mathrm{\Psi }^{}D\overline{D}`$ coupling constant is 12.84; the $`\mathrm{\Psi }^{\prime \prime }D\overline{D}`$ coupling from VMD is to 19.94, which is close to the result from direct $`\mathrm{\Psi }^{\prime \prime }D\overline{D}`$ decay given by Eq.1.

The phase space dependence of the $`\mathrm{\Psi }^{}`$ in-medium width from Eq.1 is shown by the solid line in Fig.2b). This result again can be compared with the prediction of the $`{}_{}{}^{3}P_{0}^{}`$ model given by Busetto

$`\mathrm{\Gamma }_{\mathrm{\Psi }^{}D\overline{D}}={\displaystyle \frac{\pi ^{1/2}\mathrm{\hspace{0.17em}2}^9\gamma ^2}{m_\mathrm{\Psi }\beta ^7\mathrm{\hspace{0.17em}3}^{11}}}[(q^2+m_D^2)(q^2+m_{\overline{D}}^2)]^{1/2}`$
$`\times q^3(15\beta ^22q^2)^2\mathrm{exp}(q^2/6\beta ^2)`$ (5)
and shown by the dashed line in Fig.2b). Here the calculations were done with parameters $`\beta `$ and $`\gamma `$ evaluated above. Note that at nuclear densities $`\rho 1.5\rho _0`$ both the $`{}_{}{}^{3}P_{0}^{}`$ model and VMD phase space estimates are in reasonable agreement. An exciting observation is that at normal nuclear density $`\rho _0`$=0.16 fm-3 the $`\mathrm{\Psi }^{}`$ width is $``$70 MeV, which is $``$250 times larger than that in vacuum.

Again, the $`\mathrm{\Psi }^{}`$ modification in nuclear matter can be measured through dilepton spectroscopy. As we found, the dilepton spectrum from $`\mathrm{\Psi }^{}(3686)`$ decay in matter might strongly overlap with $`\mathrm{\Psi }^{\prime \prime }(3770)`$ charmonium decay. Fig.3a) shows dilepton spectra from $`\mathrm{\Psi }^{}`$ and $`\mathrm{\Psi }^{\prime \prime }`$ decays unfolded from the production cross section. The solid lines show the results for matter at normal nuclear density, while the dashed lines indicate the spectra in free space. The measurement of a broad peak around $`\mathrm{\Psi }^{}`$ pole might be considered as an direct indication of in-medium modification of charmonia.

On the other hand, a large $`D\overline{D}`$ component of $`\mathrm{\Psi }^{}`$ charmonium should result in strong $`\mathrm{\Psi }^{}`$ absorption in nuclear matter, similar to that found Sibirtsev2  for $`J/\mathrm{\Psi }`$ dissociation. Moreover, $`\mathrm{\Psi }^{}`$ melting in nuclear matter will additionally suppress $`\mathrm{\Psi }^{}J/\mathrm{\Psi }`$ decay and partially eliminate $`J/\mathrm{\Psi }`$ yield in heavy ion collisions.

## 4 $`\chi _{c2}`$(3556).

The $`\chi _{c2}`$ charmonium level crosses the $`D\overline{D}`$ threshold at a nuclear density of about 0.28 fm-3. In vacuum the $`\chi _{c2}`$ width is 2 MeV and partial decay into the $`J/\mathrm{\Psi }`$ is $``$13.5%. There is no reliable way to estimate $`\chi _{c2}D\overline{D}`$ coupling and to evaluate the phase space dependence of the $`\chi _{c2}`$ in-medium width.

In the $`{}_{}{}^{3}P_{0}^{}`$ model the $`\chi _{c2}D\overline{D}`$ width is given as Busetto

$`\mathrm{\Gamma }_{\chi _{c2}D\overline{D}}={\displaystyle \frac{\pi ^{1/2}\mathrm{\hspace{0.17em}2}^{12}\gamma ^2}{5m_\chi \beta ^5\mathrm{\hspace{0.17em}3}^8}}[(q^2+m_D^2)(q^2+m_{\overline{D}}^2)]^{1/2}`$
$`\times q^5\mathrm{exp}(q^2/6\beta ^2),`$ (6)
where $`m_\chi `$ is the $`\chi _{c2}`$ mass and parameters $`\beta `$ and $`\gamma `$ are listed above. The in-medium $`\chi _{c2}`$ width is shown in Fig.3b) as a function of nuclear matter density. The $`\chi _{c2}`$ modification becomes significant only at large densities.

## 5 Conclusion.

Modification of the $`D`$ and $`\overline{D}`$ meson masses in nuclear matter leads to a substantial increase of the $`\chi _{c2}`$, $`\mathrm{\Psi }^{}`$ and $`\mathrm{\Psi }^{\prime \prime }`$ decay widths into the $`D\overline{D}`$ channel.

The calculations with the density independent coupling constants between the $`\mathrm{\Psi }^{}`$ and $`\mathrm{\Psi }^{\prime \prime }`$ charmonium and $`D\overline{D}`$ pair results in strong and monotonic density dependence of the $`\mathrm{\Psi }^{}`$ and $`\mathrm{\Psi }^{\prime \prime }`$ in-medium widths due to increase of the final state phase space. Within the $`{}_{}{}^{3}P_{0}^{}`$ model these couplings are also considered as a function of the in-medium $`D`$ and $`\overline{D}`$ masses and as a result the $`\mathrm{\Psi }^{}`$ and $`\mathrm{\Psi }^{\prime \prime }`$ charmonium widths do not increase monotonically with nuclear density, but saturate at $`\rho _0`$=0.16 fm-3. It was found that the saturation limits are $`\mathrm{\Gamma }_{\mathrm{\Psi }^{}D\overline{D}}`$70 MeV and $`\mathrm{\Gamma }_{\mathrm{\Psi }^{\prime \prime }D\overline{D}}`$90 MeV, which can be compared with vacuum widths of 0.277 MeV and 23.6 MeV, respectively.

The $`\chi _{c2}`$ charmonium level crosses the $`D\overline{D}`$ threshold at $`\rho `$ 1.25$`\rho _0`$ and its width increases with nuclear matter density. At $`\rho `$3$`\rho _0`$ the $`\chi _{c2}`$ decay width into $`D\overline{D}`$ is about 11 MeV, which may be compared with total $`\mathrm{\Gamma }_{\chi _{c2}}`$=2 MeV width in vacuum.

We conclude that modification of the $`D`$ and $`\overline{D}`$ in matter most dramatically affects the $`\mathrm{\Psi }^{}`$ charmonium. This can be detected by dilepton spectroscopy as the appearance of a broad peak near the $`\mathrm{\Psi }^{}`$ pole or as anomalous $`\mathrm{\Psi }^{}`$ suppression in nuclear matter.

Our results are in agreement with the most recent findings Ko  that both the mass and width of $`\mathrm{\Psi }^{\prime \prime }`$(3770) depend significantly on the $`D`$ meson mass.


text:
# The Partition Function and Level Density for Yang-Mills-Higgs Quantum Mechanics

## I Introduction

The discovery of chaoticity in the classical Yang-Mills (YM) equationsref1  has attracted attention to the system of coupled quartic oscillators with the potential $`x^2y^2`$, where $`x`$ and $`y`$ are functions of time only. This system is the simplest limiting case for the homogeneous YM equations (the so-called YM Classical Mechanics) with $`n=2`$ degrees of freedom. Despite its simplicity, the $`x^2y^2`$ model exhibits a rich versatile chaotic behaviour and belongs to the most chaotic potential systems known. Not surprisingly, this potential has been used in many fields, including chemistry, astronomy, astrophysics, and cosmology (chaotic inflation).

From the quantum mechanical point of view, this system (YM Quantum Mechanics) has a discrete spectrumref2 ; ref3  in spite of having an infinite volume of energetically accessible phase space $`\mathrm{\Gamma }`$ ref4 ; ref5

$$\mathrm{\Gamma }=𝑑x𝑑y𝑑\dot{x}𝑑\dot{y}\delta \left(\frac{1}{2}\dot{x}^2+\frac{1}{2}\dot{y}^2+\frac{g^2}{2}x^2y^2E\right),$$

where the dot stands for $`d/dt`$. This potential, therefore, violates the Weil’s lawref6 , a semiclassical relation, which states that the average number $`N(E)`$ of quantum energy levels with energy less than $`E`$ is asymptotically proportional to $`\mathrm{\Gamma }`$.

Some time ago, there appeared two important papersref7 ; ref8  devoted to the calculation of the partition function $`Z(t)`$ and the asymptotic integrated level density $`N(E)`$ for the $`x^2y^2`$ model. Their method was based on an adiabatic separation in the partition function’s dependence on $`x`$ and $`y`$ out in the narrow channels of the equipotential surface $`|xy|`$ = constant. Remarkably the dependence on the boundary dividing the two regions (central ($`|x|Q,|y|Q`$) and channel ($`Q|x|,|y|<\mathrm{}`$)) with quite different physics (essentially classical for the first region and intrinsically quantum for the second region) disappears in the final result for the partition function. This insensitivity to the dependence on the boundary (the value of $`Q`$) may seem to bode well for this method of calculation.

However, a priori we do not expect that the partition function of a non-integrable system with infinite phase space volume is calculable without any approximation. Therefore we seek an alternative approach which brings to fore the fact that $`\mathrm{\Gamma }=\mathrm{}`$ for the $`n=2`$ case. With this in mind, in this paper we calculate the partition function and the asymptotic integrated level density for the so-called Yang-Mills-Higgs Quantum Mechanics (YMHQM)ref9 ; ref4 . It is interesting to examine how $`Z(t)`$ and $`N(E)`$ behave in the limit the Higgs coupling to the $`x`$ and $`y`$amplitudes vanishes and one recovers the pure $`x^2y^2`$-system. We also consider the YMHQM for the case of $`n=3`$ for which $`\mathrm{\Gamma }`$ is finiteref4 ; ref5  even for the pure YM system.

## II Yang-Mills-Higgs Mechanics

There are several mechanisms that can suppress the classical chaoticity of chaotic systems with dimension $`n=2`$ and, in particular, to the the YM system (see Ref.ref5 ). One of them is the well-known Higgs mechanism. For spatially homogeneous fields, if only the interaction of the YM fields with the Higgs vacuum is considered, the classical Hamiltonian density for $`n=2`$ is given by

$$H=\frac{1}{2}(\dot{x}^2+\dot{y}^2)+\frac{g^2}{2}x^2y^2+\frac{v^2}{2}(x^2+y^2).$$
(1)
It is knownref9  that there is a classical “phase transition” from chaos to regular motions as the vacuum expectation value of the Higgs field $`\varphi =v`$ gets large enough. Here there is one control parameter $`\kappa =\frac{g^2v^4}{H}`$, where $`H`$ is the conserved energy density. At large $`\kappa `$ the motion is regular, whereas at $`\kappa =0`$ one deals with the developed chaos of YM Classical Mechanics. In fact, chaos appears already at $`\kappa 0.60`$.ref9  Therefore, in the study of the quantum counterpart of Eq. (1), we expect that there is a transition from one type of $`N(E)`$ to another, depending on the parameter $`v`$. The analogous transition in the adjacent energy level spacing distribution, as a function of $`v`$, was predictedref10  and established in several papersref11 .

In the next section we calculate the partition function for the quantum-mechanical counterpart of the classical Hamiltonian density Eq. (1). In section IV, we calculate the integrated level density. Section V is devoted to the $`n=3`$ case. We conclude in the last section with some discussions involving a related system with potential of the form $`\frac{g^2}{2}x^2y^2+\frac{b^2}{4}(x^4+y^4)`$.

## III Partition function for YMH quantum mechanics

The quantum Hamiltonian corresponding to Eq. (1) is given by (mass $`m=1`$ in the following)

$$=\frac{\stackrel{}{𝐩}^2}{2}+\frac{g^2}{2}x^2y^2+\frac{v^2}{2}(x^2+y^2),$$
(2)
where $`\stackrel{}{𝐩}=i\mathrm{}\stackrel{}{}`$ is the momentum operator. First a word on units. All quantities below are given in units of energy $`E`$: $`[]=E`$, $`[t]=E^1`$, $`[x]=[y]=E^{1/4}`$, $`[g]=E^0`$, $`[v]=E^{1/4}`$, $`[\mathrm{}]=E^{3/4}`$. It is obvious that the operator Eq. (2) has a discrete spectrum.

Here we calculate the trace of the heat kernel $`\mathrm{exp}(t)`$, the partition function for the Hamiltonian operator. The trace is defined for any quantum operator $`𝒜_W`$ in the Wigner representation (i.e., in the classical 2n-dimensional phase space)

$$Tr(𝒜_W)=\frac{1}{(2\pi \mathrm{})^n}𝑑\stackrel{}{𝐩}𝑑\stackrel{}{𝐪}A_W(\stackrel{}{𝐩},\stackrel{}{𝐪}).$$
(3)
Thus, we have

$$Z(t)=Tr(e^t)=_0^{\mathrm{}}𝑑Ee^{tE}\rho (E).$$
(4)
The second equation expresses the fact that the partition function $`Z(t)`$ and the density of eigenstates $`\rho (E)`$ form a Laplace transform pair. The asymptotic integrated density of states $`N(E)`$ is given by the inverse Laplace transform of $`Z(t)/t`$

$$N(E)=_0^E𝑑E^{}\rho (E^{})=L^1\left(\frac{Z(t)}{t}\right).$$
(5)
Integrating Eq. (4) over $`p_x`$ and $`p_y`$ for the Hamiltonian given by Eq. (2), we get

$$Z(t)=\frac{1}{(2\pi \mathrm{})^2}\frac{2\pi }{t}𝑑x𝑑ye^{t[v^2(x^2+y^2)+g^2x^2y^2]/2}.$$
(6)
After the integration over $`x`$, we can make a change of variable from $`y`$ to $`s=\frac{1}{2}tv^2y^2`$ to obtain

$$Z(t)=\frac{1}{(2\pi \mathrm{})^2}\left(\frac{2\pi }{t}\right)^{3/2}\sqrt{\frac{2}{tv^4}}_0^{\mathrm{}}𝑑se^ss^{1/2}\left(1+\frac{s}{2z}\right)^{1/2},$$
(7)
with $`z\frac{tv^4}{4g^2}`$. But aside from a factor of $`\sqrt{2z}e^z`$, the integral is just a representation of $`K_0(z)`$, the modified Bessel function of the third kind or the MacDonald function of order $`0`$. (see, e.g., p.140 of Ref.ref12 ). Thus we arrive at the precise expression for the partition function corresponding to the Hamiltonian Eq. (2)

$$Z(t)=\frac{1}{\sqrt{2\pi }}\frac{1}{g\mathrm{}^2t^{3/2}}\mathrm{exp}\left(\frac{tv^4}{4g^2}\right)K_0\left(\frac{tv^4}{4g^2}\right).$$
(8)
So far we have not made any calculational approximations. For two uncoupled oscillators (corresponding to $`g=0`$) Eq. (8) gives (with the aid of $`K_0(z)\sqrt{\frac{\pi }{2z}}e^z`$ as $`z\mathrm{}`$)

$$Z(t)=\frac{1}{(\mathrm{}vt)^2},$$
(9)
as it should. For the sake of completeness, let us also give the partition function for one oscillator ($`n=1`$)

$$Z(t)=\frac{1}{\mathrm{}vt}.$$
(9)
From Eq. (8), it is obvious that $`Z(t)`$ diverges logarithmically in the limit $`v0`$ (the $`x^2y^2`$ model). Indeed, using $`K_0(z)\mathrm{log}(z/2)C`$ (with $`C`$ being the Euler constant) as $`z0`$, we get, for $`v0`$,

$$Z(t)\frac{1}{\sqrt{2\pi }}\frac{1}{g\mathrm{}^2t^{3/2}}\left(\mathrm{log}\frac{8g^2}{tv^4}C\right).$$
(10)
The impossibility to disentangle the coupled oscillators from the uncoupled ones is a reflection of the logarithmic divergence of the phase space volume $`\mathrm{\Gamma }`$ on the $`E`$-shell at $`v=0`$ for $`n=2`$, as we have mentioned above. For $`n=3`$ (see section V below), $`\mathrm{\Gamma }`$ is finite and the precise $`v=0`$ limit exists.

## IV Integrated level density $`N(E)`$ ($`n=2`$)

It is well-known that the asymptotic integrated density of states $`N(E)`$ is related to the small $`t`$ divergence of the partition function $`Z(t)`$, according to the Karamata-Tauberian theorem (see, e.g., Ref. ref3 ). Using Eq. (8) and Eq. (5), we obtain the precise expression for $`N(E)`$ (with the aid of formula 3.16.3.3 of Ref.ref13 )

$`N(E)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{1}{g\mathrm{}^2}}L^1\left({\displaystyle \frac{1}{t^{5/2}}}\mathrm{exp}\left({\displaystyle \frac{tv^4}{4g^2}}\right)K_0\left({\displaystyle \frac{tv^4}{4g^2}}\right)\right)`$ (11)
$`=`$ $`{\displaystyle \frac{E^2}{2\mathrm{}^2v^2}}F({\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}};3;{\displaystyle \frac{2g^2E}{v^4}}),`$
where F is the Gauss hypergeometric function. For $`g=0`$ (the case of two uncoupled oscillators) we have, from Eq. (11),

$$N(E)=\frac{1}{2\mathrm{}^2v^2}E^2,$$
(12)
in agreement with Eq. (9). For completeness, for one oscillator, we have, from Eq. (9),

$$N(E)=\frac{1}{\mathrm{}v}E.$$
(12)
For the $`v0`$ limit, one needs the asymptotic expression of $`F(a,b;c;z)`$. Using formula 9.7.7 of Ref.ref12 , we get

$$N(E)\frac{2\sqrt{2}}{3\pi }\frac{E^{3/2}}{g\mathrm{}^2}\left(\mathrm{log}\frac{g^2E}{v^4}+5\mathrm{log}2\frac{8}{3}\right).$$
(13)
One can verify that the $`Z(t)`$ and $`N(E)`$ for Eq. (9), Eq. (9), and the corresponding Eq. (12), Eq. (12), and the ratio of the logarithmic parts of Eq. (10) and Eq. (13), are in accordance with Theorems 1.1 and 1.4 of Ref. ref3 . Also from Eq. (12) and Eq. (13), we see that the transition from order ($`g=0`$) to chaos ($`g0`$ and large, $`v`$ small, $`\kappa 0`$)333We need to emphasize that, from the results of Ref.ref9 , the chaotic regime for $`N(E)`$ actually begins at finite $`v=(\kappa E/g^2)^{1/4}`$. corresponds to the change in the $`E`$-dependence of $`N(E)`$ from $`N(E)E^2`$ to $`N(E)E^{3/2}\mathrm{log}E`$, analogous to the change in the neighbour level spacing distribution from the Poisson distribution to the Wigner-Dyson distribution. It can be seen that as the power $`\alpha `$ of the homogeneous potential $`|xy|^\alpha `$ increases from 1 to $`\mathrm{}`$, there is ref3  a systematic decrease in the power of $`E`$ in $`N(E)`$ from 2 to 1, the latter corresponding to the case of the hyperbola billiard with $`N(E)E\mathrm{log}E`$ (see Ref.ref14 ). This observation may be relevant to the Hilbert-Polya-Berry program in identifying a quantum (chaotic) Hamiltonian whose eigenvalues reproduce the Riemann zeta-function zeros.

## V Some remarks on the three dimensional YM and YMH quantum mechanics

Let us now consider the YMH Quantum Mechanics for $`n=3`$ with the Hamiltonian (with mass $`m=1`$ again)

$$=\frac{1}{2}\stackrel{}{𝐩}^2+\frac{g^2}{2}(x^2y^2+y^2z^2+z^2x^2)+\frac{v^2}{2}(x^2+y^2+z^2).$$
(14)
For the partition function, the integrations over $`\stackrel{}{𝐩}`$ can be easily carried out. Introducing the cylindrical coordinates ($`x=r\mathrm{cos}\varphi ,y=r\mathrm{sin}\varphi ,z=z`$) and performing the integrations over $`z`$ and $`\varphi `$, we obtain

$$Z(t)=\frac{1}{(2\pi \mathrm{})^3}\left(\frac{2\pi }{t}\right)^22\pi _0^{\mathrm{}}\frac{rdr}{\sqrt{v^2+g^2r^2}}\mathrm{exp}\left(\frac{tv^2r^2}{2}\frac{tg^2r^4}{16}\right)I_0\left(\frac{tg^2r^4}{16}\right),$$
(15)
where $`I_0(z)`$ is the modified Bessel function of the first kind of order $`0`$.

We see that the integral in Eq. (15) is well behaved at $`v=0`$; it is equal to

$$\frac{1}{4g}_0^{\mathrm{}}𝑑u\frac{1}{u^{3/4}}e^{g^2tu/16}I_0\left(\frac{g^2tu}{16}\right).$$
(15)
The integration over $`u`$ can be done by using formula 2.15.3.3 of Ref.ref15  to yield

$$Z(t)=\frac{\mathrm{\Gamma }^3(\frac{1}{4})}{2^{7/4}\pi ^{3/2}}\frac{1}{g^{3/2}\mathrm{}^3t^{9/4}},$$
(16)
which agrees with Eq. (4.5) of Ref.ref8  (after taking into account the extra factor of $`\frac{1}{2}`$ in the definition of the potential in Eq. (14)). The difference between the $`n=2`$ and $`n=3`$ cases is due to the fact that for $`n=3`$ the energetically accessible region pinches as $`1/x^2`$ at large $`x`$, in contrast to the $`n=2`$ case for which it pinches as $`1/x`$, leading to the logarithmic divergent result (Eq. (10)) at $`v=0`$. Note that the contibution from the channels ($`t^{1/4}x>>1`$) for $`n=3`$ is negligible in the method of separating the central part from the channels used in Ref.ref7 ; ref8 . Thus for $`n=3`$ the coupled oscillators are disentangled from the free oscillators. From Eq. (16), it follows that

$$N(E)=\frac{16}{45}\frac{2^{1/4}\mathrm{\Gamma }^2(\frac{1}{4})}{\pi ^{3/2}}\left(\frac{E^{3/4}}{\sqrt{g}\mathrm{}}\right)^3,$$
(17)
which agrees with the result given in Ref.ref8  (the first of Eq. (5.19)).

Finally, for the sake of completeness, here are the expressions for $`Z(t)`$ and $`N(E)`$ at $`g=0`$ for the case of $`n=3`$. Eq. (15) yields

$$Z(t)=\frac{1}{(\mathrm{}vt)^3}.$$
(18)
Of course, Eq. (9), Eq. (9), and Eq. (18) are in accordance with the well-known relation of the trace of the n-harmonic-oscillator heat kernel, viz., $`Z(t)=\left(2\mathrm{sinh}\frac{\mathrm{}vt}{2}\right)^n`$ in the quasiclassical approximation ($`t0`$)ref7 ; ref8 . From Eq. (18) we have, for $`N(E)`$,

$$N(E)=\frac{1}{6\mathrm{}^3v^3}E^3.$$
(19)
Eq. (16) - Eq. (19) show that the $`Z(t)`$ and $`N(E)`$ for the $`n=3`$ case are in full agreement with Theorem 1.1 of Ref.ref3 .

## VI Discussions

From the previous two sections, we see that the $`\mathrm{log}E`$ factor in $`N(E)`$ is specific to the $`x^2y^2`$-model (with or without the Higgs term), and is not generic to chaotic systems with potentials $`|xy|^\alpha `$ ($`\alpha =2,3,\mathrm{}`$) for $`n>2`$. But one may wonder whether it is generic to other homogeneous potentials for $`n=2`$. The answer depends on the type of homogeneous potentials added to the $`x^2y^2`$ term. Consider the following example of $`n=2`$ system (with finite phase space), described by the following quantum Hamiltonian

$$=\frac{\stackrel{}{𝐩}^2}{2}+\frac{g^2}{2}x^2y^2+\frac{b^2}{4}(x^4+y^4).$$
(20)
Like $`g`$, $`b`$ is dimensionless. It is known that for $`b^2=g^2/3`$ the system Eq. (20) is integrable.ref16  In the calculation of $`Z(t)`$, the integrations over $`x`$ and $`y`$ can best be done by using the polar coordinates. After the gaussian integration over the radial coordinate, the integration over the angle yields the complete elliptic integral of the first kindref12  which is related to the Gauss hypergeometric function. The result is given by

$$Z(t)=\frac{\sqrt{\pi }}{2}\frac{1}{b\mathrm{}^2t^{3/2}}F(\frac{1}{2},\frac{1}{2};1;\frac{b^2g^2}{2b^2}).$$
(21)
The corresponding asymptotic integrated density of states $`N(E)`$ is given by

$$N(E)=\frac{2E^{3/2}}{3b\mathrm{}^2}F(\frac{1}{2},\frac{1}{2};1;\frac{b^2g^2}{2b^2}).$$
(22)
For the $`b0`$ limit, we use formula 9.7.7 of Ref.ref12  again to get

$$Z(t)\frac{1}{\sqrt{2\pi }}\frac{1}{g\mathrm{}^2t^{3/2}}\mathrm{log}\frac{8g^2}{b^2}.$$
(23)
We note the logarithmic divergence as $`b0`$, i.e., when one tries to disentangle the coupled quartic oscillators from the free quartic oscillators. A simple inverse Laplace transform yields the level density

$$N(E)\frac{2\sqrt{2}}{3\pi }\frac{E^{3/2}}{g\mathrm{}^2}\mathrm{log}\frac{8g^2}{b^2}.$$
(24)
Just as there is no $`\mathrm{log}t`$ dependence in $`Z(t)`$, we find no $`\mathrm{log}E`$ dependence in $`N(E)`$ for this system. This result is not surprising from the viewpoint of dimensional considerations since both couplings in this example are dimensionless. In the case of YMHQM, the coupling $`v`$ is dimensionful, and due to this, $`\mathrm{log}t`$ and $`\mathrm{log}E`$ appear in Eq. (10) and Eq. (13) respectively in the $`v0`$ limit.

For completeness we also give the partition function for the $`g=0`$ case which corresponds to an integrable system. Using $`F(\frac{1}{2},\frac{1}{2};1;\frac{1}{2})=\mathrm{\Gamma }^2(\frac{1}{4})/(2\pi ^{3/2})`$, we have

$$Z(t)=\frac{\mathrm{\Gamma }^2(\frac{1}{4})}{4\pi }\frac{1}{b\mathrm{}^2t^{3/2}}.$$
(25)
The corresponding asymptotic integrated density of states is given by

$$N(E)=\frac{\mathrm{\Gamma }^2(\frac{1}{4})}{3\pi ^{3/2}}\frac{E^{3/2}}{b\mathrm{}^2}.$$
(26)
For another integrable case with $`b^2=g^2/3`$ in this exampleref16 , using $`F(\frac{1}{2},\frac{1}{2};1;1)=\mathrm{\Gamma }^2(\frac{1}{4})/(2\pi )^{3/2}`$, we get

$$Z(t)=\frac{\mathrm{\Gamma }^2(\frac{1}{4})}{2^{5/2}\pi }\frac{1}{b\mathrm{}^2t^{3/2}},$$
(27)
and

$$N(E)=\frac{\mathrm{\Gamma }^2(\frac{1}{4})}{3(2\pi ^3)^{1/2}}\frac{E^{3/2}}{b\mathrm{}^2}.$$
(28)
Finally, for the third integrable case with $`b=g`$ we immediately obtain the $`Z(t)`$ and $`N(E)`$ from Eq. (21) and Eq. (22) using $`F(\frac{1}{2},\frac{1}{2};1,0)=1`$.

## Acknowledgements

This work was supported in part by the US Department of Energy and by the Bahnson Fund of University of North Carolina at Chapel Hill. We thank Berndt Muller for his interest in this work, and Matthias Brack for helpful correspondence and for calling our attention to related work given in Ref. 17.