format: text

text:
# Inter-galactic Shock Acceleration and the Cosmic Gamma-ray Background

## 1 Introduction

After the re-ionization epoch, the inter-galactic medium is heated up to temperatures $`T10^4`$ K by energy input from quasars and star forming regions. As the wavelength of primordial density perturbations entering the non-linear stage, $`\lambda _{NL}(z)`$, grows longer in size, the temperature of the low red-shift gas is further raised by shock waves. These convert the specific kinetic energy associated with accretion flows, of order $`\lambda _{NL}^2(z)H^2(z)`$ \- with $`H(z)`$ the Hubble parameter and $`z`$ the cosmological red-shift - into thermal energy corresponding to 7-8 keV per baryon in massive clusters of galaxies and below $`1`$ keV in the diffuse inter-galactic medium (Cen & Ostriker, 1999). The importance of this process in cosmic history is emphasized by the fact that on the one hand the large scale structure is most prominent in thermal X-ray radiation by hot intra-cluster gas. On the other hand numerical simulations predict that most of the baryons, which are observed at high red-shift as Ly-$`\alpha `$ absorbers, end up shock-heated to temperatures in the range $`10^510^7`$ K in the contemporary universe (Cen & Ostriker, 1999).

Shock waves in cosmic environment are collision-less meaning that the dissipation process occurs via random electro-magnetic fields in the turbulent post-shock flow. In addition to a thermal population of hot gas, they generate a supra-thermal distribution of high energy particles, generally referred to as cosmic-rays (CR). Shock acceleration of CR electrons up to TeV energies is directly observed in supernova remnants (Allen et al., 1997; Koyama et al., 1997; Tanimori et al., 1998; Muraishi et al., 2000; Allen et al., 2001) and the bulk of the Galactic CR ion population up to $`10^{14}10^{15}`$ eV is also thought to originate there. The same acceleration mechanism is expected at work at large scale structure shocks. Indeed, we do observe non-thermal radiation in several massive clusters of galaxies both at radio (e.g., Feretti, 1999, and references therein) and hard X-ray frequencies (Fusco-Femiano et al., 1999, 2000, and references therein). At least in some models CR shock acceleration of both ions and electrons seems an essential ingredient to explain different features of the observed radio emission (e.g., Miniati et al., 2001a). In addition, the possibility that CRs be accelerated at inter-galactic shocks has raised the issue that CR ions, due to their long lifetime against energy losses and to their efficient confinement by magnetic irregularities (Völk et al., 1996; Colafrancesco & Blasi, 1998), may accumulate inside forming structures and store a significant fraction of the total pressure there (Miniati, 2000; Miniati et al., 2001b). Beside cosmic shocks, the energetics characterizing other CR sources such as termination shocks driven by galactic winds (Völk et al., 1996) and radio galaxies (Enßlin et al., 1997) are also noteworthy.

Along these lines of observational and theoretical progress, the idea has recently been put forth that CR electrons accelerated at inter-galactic shocks may generate the observed $`\gamma `$-ray background (CGB) by way of IC (IC) scattering off cosmic microwave background (CMB) photons (Loeb & Waxman, 2000). The CGB was measured by EGRET onboard the Compton Gamma Ray Observatory. It consists of a differential flux, $`\epsilon ^2J_\epsilon `$ 5.4 ($`\epsilon `$ /keV)-0.1 keV cm-2 s-1 sr-1 (Sreekumar et al., 1998) throughout an energy range extending from about 100 MeV to 100 GeV. The integrated photon flux above 100 MeV corresponds to 1.45 $`\times 10^5`$ ph cm-2 s-1 sr-1 (Sreekumar et al., 1998), in agreement with previously measures carried out by the SAS 2 satellite (Thompson & Fichtel, 1982). Its origin is most intriguing. The contribution from identified EGRET blazars to the integrated photon flux above 100 MeV amounts to 1$`\times 10^6`$ ph cm-2 s-1 sr-1 (Mukherjee & Chiang, 1999); that is $``$ 7 % of the measured CGB flux in corresponding units. A number of authors have suggested that the contribution from an unresolved, faint end of the blazar $`\gamma `$-ray distribution function could fill the gap between the observed flux and that which is already accounted for. And in fact, though with large uncertainties, some authors were able to validate the above conjecture by using a $`\gamma `$-ray luminosity function derived from the blazar radio luminosity function and the assumption of a constant ratio of radio to $`\gamma `$-ray fluxes (e.g., Padovani et al., 1993; Stecker & Salamon, 1996). However, the correlation between the blazar $`\gamma `$-ray and radio flux is only weakly established and the existing evidence might indeed be dominated by bias effects (Mücke et al., 1997). In fact, based solely on the observationally derived luminosity function of $`\gamma `$-ray loud AGNs that have been detected by the EGRET experiment, Chiang & Mukherjee (1998) found that only 25 % of the CGB can be accounted for by undetected $`\gamma `$-ray loud blazars. Alternatively, the idea was put forward that the CGB is generated by the decay of neutral $`\pi `$-mesons, $`\pi ^0\gamma \gamma `$, generated by an all pervasive population of CR ions interacting via p-p collisions with the nuclei of the diffuse inter-galactic medium (Dar & Shaviv, 1995) . However, it was later argued that the CR ions from the available sources of CRs in the IGM, namely shocks, Active Galactic Nuclei and normal galaxies, could provide at most a fraction of order of a % of the observed CGB (Berezinsky et al., 1997; Colafrancesco & Blasi, 1998).

In this paper we carry out a numerical study to investigate the contribution to the CGB from CRs accelerated at cosmological, inter-galactic shocks associated with structure formation. We show that the heating of the IGM is due primarily to moderately strong shocks, which allow for some fraction of the shock ram pressure to be converted into CR pressure. According to our results, the contribution to the CGB from cosmological CRs could indeed be significant. It is dominated by IC emission from primary electrons which, however, we estimate to be smaller than previously found (Loeb & Waxman, 2000). In addition we find a $`\gamma `$-ray flux from $`\pi ^0`$-decay that is not negligible and larger than, although of the same order as, the upper limits estimated by Colafrancesco & Blasi (1998).

The paper is organized as follows: we describe the numerical model in §2, present the results in §3 discuss them in §4 and summarize with §5.

## 2 The Simulation

Our investigation is based on numerical simulations of large scale structure formation which include shock acceleration, transport and energy losses/gains of CR ions and electrons. In the following sections we describe the salient features of the employed numerical techniques.

### 2.1 Cosmological Model

The formation and evolution of the large scale structure is computed by means of an Eulerian, grid based Total-Variation-Diminishing hydro+N-body code (Ryu et al., 1993). The code is capable of sharply capturing both very strong and weak shocks while being computationally relatively inexpensive. This feature is of primary importance in this kind of study since cosmic shocks have direct impact on the CR properties.

We focus on a single, currently favored $`\mathrm{\Lambda }`$CDM “concordance model” (Ostriker & Steinhardt, 1995). The assumed model is flat (e.g., Jaffe et al., 2001) with a total mass density $`\mathrm{\Omega }_m=0.3`$ and a vacuum energy density $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_m=0.7`$. The normalized Hubble constant is taken to be $`hH_0/100`$ km s-1 Mpc-1 = 0.67 (Freedman, 2000) and the baryonic mass density, $`\mathrm{\Omega }_b=0.04`$. This value, after Big Bang nucleosynthesis calculations (e.g., Olive et al., 2000) is only marginally consistent (at the 2 $`\sigma `$ level) with the best fit to cosmic microwave background (CMB) measurements (e.g., Jaffe et al., 2001). The initial density perturbations are generated as a Gaussian random field with a power spectrum characterized by a spectral index $`n_s=1`$ and “cluster-normalization” $`\sigma _8=0.9`$. The initial velocity field is computed through the Zel’dovich approximation. The computational box is set to a comoving size $`L=50`$ $`h^1`$Mpc. The dark matter component is described by 2563 particles whereas the gas component is evolved on a comoving grid of 5123 zones. Thus each numerical cell measures about 100 $`h^1`$kpc (comoving) and each dark matter particle corresponds to $`2\times 10^9`$ $`h^1`$M.

### 2.2 Cosmic-rays

The CR dynamics is computed numerically through the code COSMOCR (Miniati, 2001). The code follows simultaneously three CR populations, namely primary ions and electrons and secondary electrons/positrons ($`e^\pm `$). Primary ions and electrons are injected and accelerated at cosmic shocks according to the diffusive shock acceleration (DSA) mechanism (e.g., Blandford & Eichler, 1987). Secondary e± are produced in p-p inelastic collisions of CR ions with thermal IGM nuclei. In addition the code accounts for spatial transport and energy losses/gains undergone by each CR component. Before describing the numerical technique in some detail (§2.2.1 and §2.2.2) we address some of the physical assumptions underlying our model.

In general CRs are treated as passive quantities, meaning that at this stage their dynamical role is completely neglected both on the shock structure (test particle limit) and the gas dynamics. Recently developed kinetic models of DSA theory certainly allow for solutions in which shock acceleration is highly efficient and the CRs strongly affect the shock structure (e.g., Malkov & Drury, 2001). In these solutions, however, most of the shock associated ram pressure is converted into CR pressure which becomes the dominant dynamic component. Although it cannot be completely ruled out that such a situation occurs in the cosmic environment, both dynamical considerations (Markevitch & Vikhlinin, 1997; Horner et al., 1999; Nevalainen et al., 2000; Roussel et al., 2000; Miralda-Escude & Babul, 1995; Wu, 2000) as well as existing upper limits on $`\gamma `$-ray emission from nearby clusters (Blasi, 1999; Miniati et al., 2001b), suggest that the ICM pressure is dominated by a thermal component, with room for a still significant CR pressure component at the level of a few tens of percent. The latter is likely to affect the gas-dynamics of the intra-cluster (and intergalactic medium) and will be considered in future work.

Also, we implicitly assume the existence of a diffuse background magnetic field which is essential not only for the acceleration of the CRs but also for their subsequent confinement within cosmic structures. In support of this assumption we notice that inter-galactic magnetic fields are indeed observed in the media of clusters of galaxies at the $`\mu `$G level (Clarke et al., 2001). Observations are much more difficult for the case of filaments and therefore the evidence there is more sparse. However, the cases of rotation measure from the Coma super-clusters (Kim et al., 1989) as well as the recent detection of a radio filament (Bagchi et al., 2002) suggest that indeed sub-$`\mu `$G magnetic fields might be common in the IGM (Kronberg, 1994).

#### 2.2.1 Cosmic-ray Ions

The first step in simulating CRs consist in detecting shocks where their injection and acceleration takes place. Shocks are identified when a converging flow ($`𝐯<\mathrm{𝟎}`$) experiences a pressure jump $`\mathrm{\Delta }P/P`$ above a given threshold. For the present simulation, such threshold was set to the pressure jump across a shock with Mach number M=1.5. Successful detection of shocks requires a scheme that reproduces them as sharp transition in the numerical solution. This is especially important because relatively weak shocks, which tend to be harder to capture, cannot be neglected. In fact, they turn out to process a substantial fraction (but not all) of the kinetic energy associated to cosmic flows (Miniati et al., 2000, §3.1).

The simulated CR ions are supposed to be injected in the DSA mechanism from the thermal IGM. Here the injection process is modeled after the thermal leakage prescription (e.g., Kang & Jones, 1995). We assume that, upon shock passage, the gas thermalizes to a Maxwellian distribution defined by a post-shock temperature $`T_{\mathrm{shock}}`$. Thermal ions in the high energy end of such a distribution might be fast enough to “leak” back upstream of the shock against trapping by the plasma waves that moderate the shock itself (Ellison & Eichler, 1984; Kang & Jones, 1995; Malkov & Völk, 1995). These ions are thus injected from the thermal pool in the DSA mechanism and quickly gain energy much above the thermal average (e.g., Blandford & Eichler, 1987). As in Miniati et al. (2001b), based on cluster/group properties, we it found appropriate to follow CR ions up to momenta $`p_{max}=10^6`$GeV/c. In fact, CR ions up to energies of $`10^910^{10}`$ GeV/c can be accelerated at cosmic shocks (Kang et al., 1996), although their confinement within cosmic structures becomes difficult beyond $`10^7`$ GeV (Völk et al., 1996; Colafrancesco & Blasi, 1998). The DSA process takes place on time-scales much shorter than the simulation time-step and, from a numerical point of view, it is treated as “instantaneous”. Thus, the injected particles are redistributed to a supra-thermal population in accord with the shock properties. In the test-particle limit adopted here, this is a power law in momentum given by (e.g., Drury, 1983)

$$f(p)_{\mathrm{shock}}=f(p_{inj})_{\mathrm{Maxwell}}\left(\frac{p}{p_{inj}}\right)^q.$$
(1)
extending from $`p_{inj}`$ to $`p_{max}`$. The log-slope of the distribution is $`q\mathrm{ln}f/\mathrm{ln}p=3(\gamma _{gas}+1)/[2(1M^2)]=4/(1M^2)`$ for $`\gamma _{gas}=5/3`$ where $`M`$ is the shock Mach number and $`\gamma _{gas}`$ is the gas adiabatic index. The Maxwellian and the power-law components that make up the full distribution function join smoothly roughly at the injection momentum, $`p_{inj}`$; that is the momentum threshold above which ions are energetic enough to escape upstream the shock as mentioned above. The choice of $`p_{inj}`$ sets the normalization in the expression (1). Here, as in previous calculations, the momentum threshold is assumed as large as a few times the peak thermal value (e.g., Kang & Jones, 1995), namely

$$p_{inj}=c_1\mathrm{\hspace{0.17em}2}\sqrt{m_pk_\mathrm{B}T_{\mathrm{shock}}}$$
(2)
where $`m_p`$ is the proton mass, $`k_\mathrm{B}`$ is the Boltzmann constant. In this study we adopt a parameter value $`c_12.5`$, on the relatively low-efficiency side of canonical values. The sign, $`>`$, indicates that when for very strong shocks ($`M10`$) the produced CR pressure is larger than allowed by the test particle approximation, we additionally impose an upper limit on it by demanding that it is $`40`$% of the shock ram pressure. This choice is admittedly artificial, although it does resemble the effect of a CR back-reaction on the shock structure which suppresses the injection at the sub-shock (e.g., Berezhko & Ellison, 1999). More importantly, however, it allows us to control how the adopted parameter affects the final result.

In order to follow the further evolution of the injected CRs we divide momentum space into $`N_p`$ logarithmically equidistant intervals (momentum bins) and for each point mesh of the spatial grid, $`𝐱_𝐣`$, we define the following piecewise power-law distribution function (Jones et al., 1999; Miniati, 2001)

$$f(𝐱_i,p)=f_j(𝐱_i)p^{q_j(𝐱_i)},1<p_{j1}pp_j,$$
(3)
where $`p_j\mathrm{}`$ are the momentum bins’ extrema. Spatial propagation and energy losses/gains of the accelerated CRs are then followed by solving numerically a “kinetic” equation written in comoving coordinates and integrated over the aforementioned momentum bins, namely (Miniati, 2001)

$$\frac{n(𝐱_𝐣,p_i)}{t}=\left[𝐮n(𝐱_𝐣,p_i)\right]+\left[b(p)\mathrm{\hspace{0.33em}4}\pi p^2f(𝐱_𝐣,p)\right]_{p_{i1}}^{p_i}+Q(𝐱_𝐣,p_i).$$
(4)
Here $`n(𝐱,p_i)=_{p_{i1}}^{p_i}4\pi p^2f(𝐱_𝐣,p)𝑑p`$ is the number density of CR in the $`i_{th}`$ momentum bin, $`b(p)\left(\frac{dp}{d\tau }\right)_{tot}`$ describes adiabatic energy losses/gains as well as Coulomb collisions and p-p inelastic scattering and $`Q_i`$ is a source term, $`i(𝐱_𝐣,p_i)`$, integrated over that momentum bin. At each time-step, $`q_j(𝐱_i)`$ is determined self-consistently from $`n(𝐱_i,p_j)`$ and the required continuity of $`f(𝐱_i,p_j)`$ at bins interfaces (Jones et al., 1999). Here we employed 4 momentum bins, although our tests indicate that the results are unchanged when $`N_p`$ is doubled.

#### 2.2.2 Semi-implicit COSMOCR

Injection of thermal electrons in the acceleration mechanism is here simply modeled by assuming that the ratio between CR electrons and ions at relativistic energies is $`R_{e/p}=10^2`$. Observationally, for Galactic CRs this ratio has been measured in the range $`1\times 10^25\times 10^2`$ (Müller & Tang, 1987; Müller & et al., 1995; Allen et al., 2001). The introduction of this parameter (see also Ellison et al., 2000) simplifies the treatment of this process, which is otherwise quite complex. In fact, the shock dynamics is presumably dominated by the proton component with which most of the mass, momentum and energy are associated. However, for an effective electron-wave interaction the resonant condition requires $`\omega k_{}v_{}=\pm \omega _{ci}`$, where $`\omega `$ and $`k`$ are the wave frequency and wave-vector respectively, $`v`$ is the particle speed, the symbol $``$ indicates the component parallel to the ordered magnetic field, $`B_0`$, and $`\omega _{Le,i}=eB_0/m_{e,i}c`$ is the ionic/electronic ($`i,e`$) Larmor frequency. Thus, wave-modes generated by leaking protons, with $`k\omega _{Li}/v`$, are much lower than those with which a thermal electron can interact. Therefore, either electrons are pre-accelerated (e.g., McClemnts et al., 1997) or a different injection model is required (e.g., Levinson, 1992, 1994; Bykov & Uvarov, 1999).

The dynamics of CR electrons is further complicated by the severe energy losses they suffer. For this reason the numerical treatment is far more delicate than for the CR ions as the distribution function steepens and cuts off above a certain momentum. In order to determine the evolution of this component more accurately, in addition to the electron bin number density described by eq. (4), we also follow the associated kinetic energy $`\epsilon _i=4\pi _{p_{i1}}^{p_i}p^2f(p)T(p)𝑑p`$. Here $`T(p)=(\gamma 1)m_ec^2`$ is the particle kinetic energy and $`\gamma =1/\sqrt{1(v/c)^2}`$ is the Lorentz factor. The equation describing the evolution of $`\epsilon (𝐱_𝐣,p_i)`$ is obtained in analogy to eq. (4), by integrating over the $`i_{th}`$ momentum bin a kinetic equation that has been multiplied by $`T(p)`$. This reads (Miniati, 2001)

$$\frac{\epsilon (𝐱_𝐣,p_i)}{t}=\left[𝐮\epsilon (𝐱_𝐣,p_i)\right]+\left[b(p)\mathrm{\hspace{0.33em}4}\pi p^2f(𝐱_𝐣,p)T(p)\right]_{p_{i1}}^{p_i}_{p_{i1}}^{p_i}b(p)\frac{4\pi cp^3f(𝐱_𝐣,p)}{\sqrt{m_ec^2+p^2}}𝑑p+S(𝐱_𝐣,p_i),$$
(5)
where $`S(𝐱_𝐣,p_i)=4\pi _{p_{i1}}^{p_i}p^2T(p)i(𝐱_𝐣,p_i)`$ and the third term on the right hand side, so written for numerical reasons, includes a combination of contributions from CR pressure work and sink terms. Thus the evolution of the CR electrons is described by eq. (4) and (5) and now the slope of the distribution function at each grid point and momentum bin is determined self consistently by the values of $`n(𝐱_𝐣,p_i)`$ and $`\epsilon (𝐱_𝐣,p_i)`$. In addition, at each grid zone we also follow explicitly the upper momentum cut-off in the electron distribution function.

For an accurate time integration of the system of eq. (4) and (5) the explicit formulation described in Miniati (2001) required a time-step, $`\mathrm{\Delta }t`$, much shorter than the electron cooling time, $`\tau _{cool}`$. And in fact, in order to properly compute in a cosmological simulation (where $`\mathrm{\Delta }t10^710^8`$yr) the electrons responsible for synchrotron and inverse Compton emission at the frequencies of interest, Miniati et al. (2001a) had to appropriately subcycle the integration of eq. (4) and (5) for several iterations. This procedure is not computationally convenient particularly at high red-shift. Furthermore, it becomes impractical when, in order to compute cosmic background radiation as in the present instance, the CR electron population is needed throughout most of the simulation. For this reason, for the electron population we have reformulated the momentum space portion (i.e., without the advection term) of the numerical scheme in a semi-implicit form. Such reformulation is quite straightforward an exercise and it mainly consists of expressing all the contribution to the fluxes in momentum space in terms of $`n(𝐱_𝐣,p_i)`$ and $`\epsilon (𝐱_𝐣,p_i)`$ at times $`t`$ and $`t+\mathrm{\Delta }t`$. The semi-implicit formulation is nominally less precise than the explicit one, although this will only affect the region near the momentum cut-off. Here, however, even though we explicitly follow the upper momentum cut-off, the major limitations of the numerical model are due to the sub-grid power-law description of the distribution function rather than to the semi-implicit formulation. Therefore, there is no significant loss of accuracy with our new approach.

Given the finite simulation time-step, the time evolution of the very high and very low energy ends of the distribution function, characterized by a cooling time $`\tau _{cool}(p)<\mathrm{\Delta }t`$, cannot be followed even with a semi-implicit scheme. We notice that in the transport equation the diffusion term can be neglected away from shocks (Miniati, 2001; Jones et al., 1999) and CRs are basically transported by the fluid (to which they are coupled through the magnetic field) at the flow speed $`v`$. Thus, because of the Courant condition, $`v\mathrm{\Delta }t\mathrm{\Delta }x`$, where $`\mathrm{\Delta }x`$ is the mesh size, electrons such that $`\tau _{cool}(p)<\mathrm{\Delta }t`$, will not propagate outside the grid zone in which they have been created. Therefore, in these energy ranges it is appropriate to take the steady state solution to eq. (4) as

$$n(𝐱_𝐣,p_i)=4\pi _{p_{i1}}^{p_i}p^2f_c(𝐱_𝐣,p_i)𝑑p=4\pi _{p_{i1}}^{p_i}\frac{dp}{b(p)}_p^{\mathrm{}}\rho ^2i(𝐱_𝐣,\rho )𝑑\rho .$$
(6)
We caution that for the case of shock acceleration, the above steady-state solution expresses the balance between a source and a loss term within the volume of a grid zone, but not at the acceleration region. In fact, the rate of energy gain at a shock must be higher than the loss rate up to the maximum energy of the accelerated particles. Physically, the expression in eq. (6) represents a cumulative population $`f_c`$, i.e., a summation of all the individual populations of CR electrons within the grid zone $`𝐱_𝐣`$ that, after they emerge from the acceleration region, and as they are being advected away from the shock, start aging. Aging here refers to the modifications produced by energy losses. It is easy to infer from eq. (6) that if the source term is a power law with slope $`q_s`$, the distribution $`f_c`$ will have an energy dependent slope

$$q_c(p)=\frac{\mathrm{ln}f_c}{\mathrm{ln}p}=q_s1+\frac{\mathrm{ln}b(p)}{\mathrm{ln}p}.$$
(7)
For dominant synchrotron and inverse Compton losses, as in the case of very high energy electrons, eq. (7) gives $`q_cq_s+1`$, whereas when Coulomb collisions dominate $`q_cq_s1`$,. Thus, effectively because of cooling and of a finite computational time-step, the numerical solution at a shock will be a broken power-law. At high energy the break will occur at momentum

$$p_{break}^h10^4\left(\frac{\mathrm{\Delta }t}{10^8\mathrm{yr}}\right)^1\left(1+\frac{\mathrm{U}_\mathrm{B}}{\mathrm{U}_{\mathrm{CMB}}}\right)^1\mathrm{MeV},$$
(8)
where $`\mathrm{U}_\mathrm{B}`$ and $`\mathrm{U}_{\mathrm{CMB}}`$ refer to the magnetic and CMB energy density respectively. At low energy, the break corresponds to

$$p_{break}^l1.8\left(\frac{\mathrm{\Delta }t}{10^8\mathrm{yr}}\right)\left(\frac{n_{\mathrm{gas}}}{10^3\mathrm{cm}^3}\right)\mathrm{MeV}.$$
(9)
For the primary CR electrons we use 5 logarithmically equidistant bins, spanning momentum space from $`p_{min}=15`$ MeV up to $`p_2=10^2`$ GeV. As for the case of CR ions, our tests show that the final results are not affected in any significant way by increasing the number of momentum bins. In addition, the distribution function above $`p_2`$, computed according to the “steady state” solution \[eq. (6)\] for zones with non-null source terms, is recorded in an extra momentum bin ranging from $`p_2`$ up to $`2\times 10`$ TeV. The latter is an appropriate value for the maximum energy of accelerated CR electrons, provided a magnetic field of order $`0.1\mu `$G throughout the IGM (Loeb & Waxman, 2000).

Although not needed in the present study, we also have the capability of computing the sub-relativistic component below $`p_{break}^l`$. This component might be of interest, for example in order to estimate the emission of high energy radiation through non-thermal bremsstrahlung (e.g., Miniati, 2002). Notice that in our model the energy associated to both CR ions and electrons is always fraction of the shock ram pressure. The small fraction of energy that is carried by the sub-relativistic components of both these CR populations, with time is transferred back to the gas in thermal form through Coulomb collisions. However, this same energy would have been converted into thermal energy through dissipation at the shock anyway, had DSA not be present. Therefore, it is obvious that the sub-relativistic populations created in our model do not produce any “thermodynamic catastrophe” (cf. Petrosian, 2001, and references therein).

#### 2.2.3 Secondary Electrons and Positrons

At each time-step we also compute the distribution of high energy $`\gamma `$-ray emitting secondary electrons and positrons (e±) produced in hadronic collisions of CR ions with the nuclei of the intra-cluster gas, within the same momentum range as for the shock accelerated electrons. Secondary $`e^\pm `$ are generated in the decay of charged muons according to

$`\mu ^\pm e^\pm +\nu _e(\overline{\nu }_e)+\overline{\nu }_\mu (\nu _\mu )`$ (10)
which in turn are produced in the following reactions

$`p+p`$ $`\pi ^\pm +X,`$ $`\pi ^\pm \mu ^\pm +\nu (\overline{\nu })`$ (11)
$`p+p`$ $`K^\pm +X,`$ $`K^\pm \mu ^\pm +\nu (\overline{\nu })`$ (13)
$`K^\pm \pi ^0+\pi ^\pm `$
In addition to $`p+p`$ inelastic collisions, the above cascades are also triggered by the following interactions between CR ions and thermal nuclei: $`p+`$He, $`\alpha +`$H and $`\alpha +`$He. We include them by assuming a helium number fraction of 7.3% for the background gas and a ratio $`(H/He)15`$ at fixed energy-per-nucleon for the CRs (Meyer et al., 1997). Further details on the source term for secondary $`e^\pm `$ are given in Miniati (2001).

## 3 Results

The shock strength, which we classify according to the shock Mach number, $`M`$, affects the CR produced background radiation in two respects: the normalization, which turns out higher for a flatter distribution function because it allows for a larger number of high energy, $`\gamma `$-ray emitting CRs; and the spectral index which, as shown below, has a simple dependence on the distribution function log-slope $`q`$. In addition, we notice that for a given injection rate, high Mach number shocks are more efficient in converting the shock ram pressure into CR ion pressure. As a first step, therefore, we investigate the strength of the shocks that were chiefly responsible for the heating of the inter-galactic medium and from which the primary CRs originate.

The Mach number is given by the ratio of the shock speed to the local sound speed (Landau & Lifshitz, 1987). When computing the shock Mach number we assume that the pre-shock gas temperature is not less than $`10^4`$ K, the value it would quickly reach had photo-heating due to UV background light as well as radiation from the hot post-shock region been included. For this reason, the highest shock Mach numbers reported below do not exceed a few 100, contrary to the case in Miniati et al. (2000) where such lower limit on the gas temperature had not been assumed.

### 3.1 Inter-galactic Shocks

The amount of thermal energy produced at shock fronts throughout cosmic history is given by the expression

$$\mathrm{\Delta }E_{th}=_{\tau _i}^{\tau (z=0)}𝑑t𝑑S_{shock}k_B\mathrm{\Delta }T(nv)_{shock}$$
(14)
where $`t`$ is cosmic time, $`dS_{shock}`$ is an element of shock surface and $`(nv)_{shock}`$ and $`\mathrm{\Delta }T`$ the gas flux and temperature jump across it. Thus, the quantity $`dS_{shock}k_B\mathrm{\Delta }T(nv)_{shock}`$ represents the amount of kinetic energy dissipated at the shock into thermal energy per particle and per unit time. In order to determine which shocks were primarily responsible for heating the IGM, in Fig. 1 we plot the differential quantity $`\frac{^2\mathrm{\Delta }E_{th}}{\mathrm{log}M\mathrm{log}T_1}`$ as a function of pre-shock temperature, $`T_1`$, and shock Mach number, $`M`$. It shows that in terms of parameter space $`(M,T)`$, shock heating is confined to a narrow strip whose upper and lower edges are loci of constant velocity of order a few $`10^3\mathrm{km}\mathrm{s}^1`$ and $`10^2\mathrm{km}\mathrm{s}^1`$ respectively. A region where $`3<M<10`$ and $`10^6<T<\mathrm{a}\mathrm{few}\times 10^7`$ further stands out where shock heating seems most “influential”. More specifically, from a quantitative analysis we find that thermal energy due to shock heating is contributed in the amount of 30% by weak shocks ($`M<4`$), 45% by moderately strong shocks ($`4M10`$) and 25% by strong shocks ($`M>10`$). Thus we expect most of the CRs accelerated at inter-galactic shocks to be described by flat distribution functions.

### 3.2 The $`\gamma `$-ray Background

The $`\gamma `$-ray flux in units of ‘keV cm-2 s-1 sr-1’ at a given photon energy, $`\epsilon `$, is computed as

$$\epsilon ^2J(\epsilon )=\epsilon \frac{c}{4\pi H_0}_0^{z_{max}}\frac{e^{\tau _{\gamma \gamma }}}{[\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }]^{1/2}}\frac{j[\epsilon (1+z),z]}{(1+z)^4}𝑑z$$
(15)
where $`j(\epsilon ,z)`$ is the computational-box-averaged spectral emissivity in units ‘photons cm-3 s-1’ computed at each simulation red-shift, $`z`$, and at the appropriately blue-shifted photon energy $`\epsilon (1+z)`$. In addition, $`\tau _{\gamma \gamma }`$ is an attenuation factor due to photo-pair creation, $`\gamma \gamma e^\pm `$. It can actually be ignored given that the energy, $`\epsilon `$, of the measured $`\gamma `$-ray photons is below 100 GeV and that the contribution from the high red-shift universe is not significant (see below). Finally, $`c`$ indicates the speed of light and $`z_{max}`$ an upper limit of integration. The actual value of the latter parameter is not too important as most of the $`\gamma `$-ray emission is produced in the low red-shift universe, with $`j(\epsilon ,z)`$ decreasing by about an order of magnitude per unit red-shift increase. The strong red-shift evolution of $`j(\epsilon ,z)`$ is related to the available amount of kinetic energy to be dissipated at shocks, which is itself a strongly decreasing function of red-shift (Miniati et al., 2000). We consider two emission processes: IC of CR electrons scattering off CMB photons (leptonic), and decay of neutral pions, $`\pi ^0\gamma \gamma `$ (hadronic), produced in p-p inelastic collisions of CR ions off nuclei of inter-galactic gas. In addition we also consider IC emission from secondary e± produced in the aforementioned p-p collisions. Our treatment of IC emission and $`\pi ^0`$ production are detailed in Miniati et al. (2001a) and Miniati et al. (2001b) respectively. For IC emission we simply use a Thompson cross section, as Klein-Nishina corrections become relevant only for electron energies about an order of magnitude above the maximum considered here.

The results are shown in Fig. 2 where we plot $`\epsilon ^2J(\epsilon )`$ as a function of $`\epsilon `$. The total flux (solid line) corresponds to a constant value at the level of 0.2 keV cm-2 s-1 sr-1 throughout the spectrum. It is clearly dominated by IC emission from shock accelerated CR electrons (dotted line). A fraction of order 30% is also produced by $`\pi ^0`$-decay (dashed line), whereas IC emission from secondary e± (dotted-dashed line) turns out about an order of magnitude smaller, although it is of the same order as $`\pi ^0`$-decay and becomes larger than it around 100 MeV. For comparison the observational data, taken from the EGRET experiment (Sreekumar et al., 1998), are also plotted with the associated error-bars (solid dots).

Remarkably, all three components produce the same flat spectrum, similar in shape to the observed one, except below $``$ GeV where the $`\pi ^0`$ production cross section drops due to phase-space constraints. This result is a reflection of the fact that the CRs at the origin of the $`\gamma `$-ray radiation were produced in strong shocks and have a flat distribution function with a log-slope $`q4`$. Indeed, IC emissivity from a power-law distribution of CR electrons, $`f(p)p^{q_e}`$, is of the form $`j(\epsilon )\epsilon ^{(q_e3)/2}`$. Therefore, for a steady state distribution function of CR electrons accelerated at strong shocks ($`M>4`$) and subject to IC losses $`q_e5`$ \[see eq. (7)\], $`\epsilon ^2J(\epsilon )\epsilon j(\epsilon )\epsilon ^0`$. As for the hadronic component, the pion emissivity (in units of pions s-1 cm-3) from a CR ion distribution $`f(p)p^{q_p}`$, is $`j_{\pi ^0}(\epsilon _{\pi ^0})\epsilon _{\pi ^0}^{(4q_p13)/3}`$ (cf. Mannheim & Schlickeiser, 1994). Thus, for $`q_p4`$ (strong shocks), as before $`\epsilon ^2J(\epsilon )\epsilon ^0`$. Analogously, the secondary e± injected in the above hadronic processes reach a steady state distribution against IC losses of the form $`f_{e^\pm }(\epsilon _{e^\pm })\epsilon _{e^\pm }^{(4q_p1)/3}\epsilon _{e^\pm }^5`$ for $`q_p4`$ \[cf. eq. (7) and Mannheim & Schlickeiser 1994\]. As the steady state distribution of primary $`e^{}`$, this component also produces an IC emission $`j(\epsilon )\epsilon ^1`$ and a $`\gamma `$-ray background flux $`\epsilon j(\epsilon )\epsilon ^0`$.

Although with a spectral shape that accords with observations, the computed $`\gamma `$-ray flux is only 15 % of the observed CGB. A number of factors, which we address below, might affect these results. Firstly, resolution effects could be important, particularly for the hadronic component, because neutral $`\pi `$-mesons are produced in a two body process. In fact the emission from smaller structures is likely to be underestimated, although we do not expect the hadronic component to be dominated by the contribution of the low end of the mass function. This is corroborated by the fact that our present estimate for the flux produced by this process is already larger than, although of the same order of magnitude as, the upper limit given by the analytical calculations carried out by Colafrancesco & Blasi (1998). In particular, between 1 GeV and 100 GeV the computed $`\gamma `$-ray flux from $`\pi ^0`$-decay corresponds to about 5 % of the measured value and to $``$ 7% when combined with the contribution of IC emission from $`e^\pm `$. In comparison Colafrancesco & Blasi (1998) estimated a range of values $`0.52`$ % (but they did not include IC emission from secondary $`e^\pm `$). Thus the error on the $`\gamma `$-ray flux as a result of hadronic processes should amount to a factor two or so. As for shock accelerated electrons, our tests indicate that coarse grid effects are not important and that they amount at most to a few tens of percent.

More important is the issue related to the assumed acceleration efficiency, $`\theta `$, for converting the shock ram pressure into CR pressure. In this calculations we adopted $`\theta _{max}40\%`$ with only a small fraction $`R_{e/p}10^2`$ of the generated CR pressure to be borne by electrons. Shock acceleration can be even more efficient than allowed here. However, in those cases we expect non-linear effects to suppress the injection rate with consequent reduction of the $`\gamma `$-ray emitting CR population whilst most of the CR pressure is carried by the upper end of a concave-up distribution of CR ions. In principle the parameter $`R_{e/p}`$ could be higher, allowing for a substantially larger contribution from the leptonic component. This point is particularly important because it has been recently claimed that IC $`\gamma `$-ray flux from CR accelerated at inter-galactic shocks can make about all of the unexplained CGB (Loeb & Waxman, 2000). Therefore, in the following we will try to assess the efficiency of conversion of shock ram pressure into CR electrons based on EGRET observations of nearby clusters.

### 3.3 IC $`\gamma `$-ray Emission from Clusters of Galaxies

For this purpose, we have computed the integrated photon $`\gamma `$-ray luminosity above 100 MeV, L$`{}_{\gamma }{}^{}(>100`$ MeV), produced by primary CR electrons accelerated at accretion shocks associated to groups/clusters of galaxies. Lγ was computed within a volume of radius about 5 $`h^1`$Mpc around the cluster center, although most of the flux is produced roughly within a radius of 3 $`h^1`$Mpc (Miniati, 2002). Collapsed objects were identified through a slightly modified version of the “spherical over-density” method (Lacey & Cole, 1994), fully detailed in Miniati et al. (2000). However, because of the large spatial extension of the IC emission, there is a real risk of mis-association with any selected object of radiation emission produced in nearby but unrelated accretion shocks. To circumvent this problem we further refined our sample by rejecting the smaller of two or more objects that happen to be closer to each other than a distance of 7.5 Mpc.

In fig. 3 we plot the computed $`\gamma `$-ray photon luminosity, L$`{}_{\gamma }{}^{}(>100\mathrm{MeV})`$, against Tx, that is the thermal X-ray luminosity weighted cluster core temperature. The correlation between the plotted quantities is tight and the numerical data are well fit with a single power-law function. After performing a simple least-$`\chi ^2`$ analysis we obtain the best-fit parameters as

$$\mathrm{L}_\gamma (>100\mathrm{MeV})=8.7\times 10^{43}\left(\frac{R_{e/p}}{10^2}\right)\left(\frac{\mathrm{T}}{\mathrm{keV}}\right)^{2.6}\mathrm{ph}\mathrm{s}^1.$$
(16)
The temperature scaling for Lγ obtained is close to what expected based on clusters scaling relations. In fact, in a steady state configuration appropriate here due to the short CR electrons cooling time, the radiation flux must be produced at the same rate at which kinetic energy is being supplied to the external accretion shocks. That implies $`F_\gamma \rho v_{shock}^3R^2T_x^{2.5}`$. The plot in fig. 3 is noteworthy by itself. In fact it shows that the $`\gamma `$-ray flux associated to IC emission from shock accelerated electrons is comparable with the flux from $`\pi ^0`$-decay which is usually regarded as the main source of $`\gamma `$-rays from clusters. Thus, for a correct interpretation of the $`\gamma `$-ray emission from clusters that planned $`\gamma `$-ray facilities aim at detecting, this component needs to be properly addressed. This task is being undertaken in a separate paper (Miniati, 2002).

In order to constrain the electron acceleration efficiency in fig. 3 we compare the computed $`\gamma `$-ray photon luminosity above 100 MeV to upper limits set by the EGRET experiment on the $`\gamma `$-ray flux from nearby clusters of galaxies. Of the 58 cases that appear in the list by Reimer (1999), we found that the tightest constraints were those set by the experimental upper limits on the following clusters: Virgo, Perseus and Coma. The latter two have temperatures larger than those of objects found in the computational box and, therefore, their upper limits must be compared to values extrapolated from the the scaling relation given in eq. (16). The flux upper limits were converted into luminosity limits based on the relation L$`{}_{\gamma }{}^{}=4\pi d_L^2`$Fγ, where $`d_L`$ is the luminosity distance taken, for each cluster, from the literature. In table 1 we summarize the assumed values for the cluster temperature, distance and $`\gamma `$-ray flux upper limits. The observational data adopted here, although numerically very close to those published in Reimer (1999), were recently re-obtained by Reimer after careful and improved analysis of numerous staked observations for each source (Reimer, priv. comm.).

The plotted experimental upper limits (particularly those relative to Coma and Virgo clusters) require that the fraction of shock energy that goes into relativistic CR electrons be $``$0.8% $`\times [\mathrm{log}(p_{max}/m_ec)/\mathrm{log}(4\times 10^7)]`$ where $`p_{max}`$ is the maximum energy of the accelerated electrons. For the adopted ion injection efficiency this roughly translates into $`R_{e/p}2\times 10^2`$. In any case, this allows us to set an upper limit on the computed $`\gamma `$-ray flux of about 0.35 keV cm-2 s-1 sr-1. Thus, according to our computation, cosmological CRs could overall contribute an important fraction of order $`25`$ % of the extra-galactic cosmic $`\gamma `$-ray background as measured by EGRET (Sreekumar et al., 1998). Due to the scatter in the plot of fig. 3 \- which we regard as “real” and attribute to the different shock histories of the simulated objects - it is plausible to relax the upper limit on the acceleration efficiency of CR electrons by another factor $`<2`$. This is allowed, however, only if the $`\gamma `$-ray flux produced by $`\pi ^0`$-decay is negligible. And although this is an extreme scenario, it cannot be excluded because $`\gamma `$-rays from galaxy clusters/groups have never been measured. Thus in this latter case the IC $`\gamma `$-ray flux could account for up to $`30`$ % of the measured CGB. Notice that the upper limits derived here rely mainly on the assumption that we compute correctly the relative fraction of $`\gamma `$-ray flux from accretion shocks onto clusters and on smaller structures, where most of the emission actually originates.

To summarize, the whole unexplained $`\gamma `$-ray background flux could be accounted for by assuming an acceleration efficiency for the electronic component $`3\%`$, or higher by a factor $`7.5`$ than assumed in fig. 2, but that would result in fluxes from nearby clusters such as Virgo, Coma and Perseus well above the experimental upper limits. It is instructive to compare in some more detail our approach and that of Loeb & Waxman (2000) as the differences go beyond the assumed acceleration efficiency for the electronic component. In fact, in Loeb & Waxman model, the population of emitting CR electrons is computed by postulating that a fraction $`\xi _e`$ of the thermal energy in the shocked IGM is converted into a flat power-law of relativistic electrons. The thermal energy density of the shocked IGM is estimated as $`1.5f_{sh}n_pk_BT`$, where $`f_{sh}`$ is the fraction of shocked baryons, $`n_p`$ the average proton number density and $`T`$ the mass-weighted average gas temperature. For their quantitative estimates Loeb & Waxman (2000) assume $`f_{sh}(k_BT/1\mathrm{k}\mathrm{e}\mathrm{V})1`$ and $`\xi _e0.05`$. Notice that, unlike $`\xi _e`$, the efficiency parameter defined in our case is based on the shock ram pressure instead of the post-shock thermal energy density. For moderately strong and strong shocks, the two definitions differs by 12.5 % (it gets more complicated for weak shocks but this case is not relevant because weak shocks are inefficient accelerators). In particular, assuming $`\xi _e0.05`$ as in Loeb & Waxman (2000) means that 5.6 % of the ram pressure is converted into CR electrons. This corresponds to a factor $`1519`$ larger than the value of 0.3-0.4 % adopted to compute the dotted line in fig. 2 which lies, however, only a factor 10 below Loeb & Waxman ’s estimate. This difference by a factor $`1.52`$ grows if one accounts for the following additional matters: effectively in the simulation $`f_{sh}2/3`$ and $`k_BT/keV1/2`$ (see also Keshet et al., 2002). In addition, as pointed out in §3.1, only $`2/3`$ of the IGM thermal energy comes from shocks strong enough (i.e., $`M4`$) to produce $`\gamma `$-ray emitting CR electrons. Thus, using the appropriate fudge factors Loeb & Waxman expression would produce a result smaller than the simulation by almost an order of magnitude (factor $`7`$). In part this is due to the fact that, by following the full evolution of the shock accelerated CR electrons from 10 MeV to 2 TeV, we are also able to account for a re-acceleration process. That is the re-energization of aged (by cooling) CR electrons that were produced in earlier shock events. As already pointed out, this effect increases our computed CGB by a factor between 2 and 3. Part of the remaining discrepancy could be due to $`\gamma `$-ray flux emitted at high red-shift ($`z>0.5`$) that Loeb & Waxman neglected, although according to our simulation that should not be too significant. Most likely, however, it reflects the fact that the actual conversion of shock energy into CR electrons cannot be simply estimated through the IGM average temperature. In any case, given its simple formulation that estimate is quite close to our results.

Just prior to submission of this paper, the work by Keshet et al. (2002) was brought to our attention. The authors further investigate the Loeb & Waxman (2000) model utilizing smoothed particle hydrodynamic (SPH) numerical simulations. They detect the occurrence of shocks in individual SPH particles by measuring changes in their entropy every some time interval. When shocks are identified a fraction $`\xi _e`$ of the generated post-shock thermal energy is converted into CR electron energy. The CR electrons are distributed in energy according to the test particle limit of DSA theory, as assumed here. Their results, particularly the level of CGB produced by shock accelerated electrons, is in good agreement with our estimated flux plotted in fig. 2. However, there are some differences that are worth pointing out. In fact, for the same estimated $`\gamma `$-ray flux our simulation requires an electron acceleration efficiency smaller than theirs by about an order of magnitude (of course this includes already the different definitions of efficiency). We find that most of the thermal energy in the low red-shift universe was produced by what we refer to as moderately strong shocks, that is shocks with Mach number in the range 4-10. Those authors mention some difficulty in tracking these shocks. If these are neglected according to our results in §3.1 one would be considering only about 1/3 of the total flux through shocks (although their estimate indicates more of a factor 1/2). In addition our simulation also includes a re-acceleration process which, as mentioned before, boosts the computed flux by a factor 2-3. These two factors, namely that related to moderately strong shocks and re-acceleration, are probably the most important in determining the differences between our and Keshet et al. results.

### 3.4 Emission Isotropy

In order to assess the spatial distribution of the $`\gamma `$-ray emissivity, in fig. 4 we show a histogram of the total integrated photon flux above 100 MeV, F$`{}_{\gamma }{}^{}(>100`$ MeV)= $`_{100\mathrm{M}\mathrm{e}\mathrm{V}}J(\epsilon )𝑑\epsilon `$, versus the temperature of the inter-galactic medium where it originates for the three computed components. The diagram illustrates that the contributions from both $`\pi ^0`$-decay (solid line) and secondary e± (dashed line) are strongly concentrated within regions with temperature $`T10^7`$keV, corresponding to groups and clusters of galaxies. On the contrary the IC emission by shock accelerated CR electrons (dotted line) is roughly equally produced in the shocked IGM with temperatures ranging from a few $`\times 10^5`$ K up to several $`\times 10^7`$ K. Notice, however, that because of strong radiative losses, these high energy electrons are almost exclusively found in “narrow” post-shock regions. Thus for this component (dotted curve), the temperature axis in fig. 4 is indicative of the gas temperature behind shocks, rather than anywhere throughout the diffuse IGM. In any case, the radiation flux produced by this component is expected to be highly isotropic (Waxman & Loeb, 2000) and, as it originates from the shocked diffuse IGM, if detected will allow an independent lower limit on the cosmic baryon density (Loeb & Waxman, 2000).

For all components, the $`\gamma `$-ray surface brightness will turn out most prominent toward galaxy clusters even for the IC emission. Waxman & Loeb (2000) predict for the IC component fluctuations of the $`\gamma `$-ray background intensity of order 40 % on the sub-degree scale. This fluctuations should actually be even larger when the hadronic component that they neglected is accounted for. In fact, the latter should be apparent only toward massive clusters and should dominate the emission within the cluster core (Miniati, 2002).

## 4 Discussion

We have constrained the energy acceleration efficiency of relativistic electrons to $`1\times [\mathrm{log}(p_{max}/m_ec)/\mathrm{log}(4\times 10^7)]\%`$ based on EGRET experimental upper limits on the photon flux above 100 MeV from nearby clusters. We specify that the found efficiency limit is valid only insofar as the accelerated CR electrons are power-law distributions extending above $`p_{max}\sqrt{75\mathrm{M}\mathrm{e}\mathrm{V}/\epsilon _{\mathrm{CMB}}}`$ mec $`100`$ GeV, where $`\epsilon _{\mathrm{CMB}}`$ is the average CMB photon energy. And, in fact, strictly speaking the found constraint actually regards the total number of CR electrons in this energy range. This upper limit on the electron acceleration efficiency plays a crucial role in the determination of the allowed contribution from cosmological CRs to the CGB. One could doubt that the EGRET measurements would apply for extended sources such as galaxy clusters accretion shocks, which measure several degree on the sky. However, the EGRET field of view is much larger than this and the exposure times drop by only a factor of two at about 18$`\mathrm{°}`$ from the instrument axis. Thus, had there been any extended flux above the upper limit set by the EGRET experiment they would have been detected (O. Reimer priv. comm.).

The upper limit just obtained on electron acceleration efficiency strictly speaking concerns the acceleration at cluster shocks, whereas a significant part of the $`\gamma `$-ray emission is produced in inter-galactic shocks. However, since clusters of galaxies are the deepest potential wells in cosmic environment, their associated shocks are indeed the strongest. Therefore, if anything, the emissivity is expected to be highest there, unless these shocks are strongly modified by a back-reaction of CR pressure. In addition, unlike for massive clusters, shocks from accretion flows on filaments and smaller structures in general could be weakened by energy input into the IGM from supernovae, stellar winds and radio-galaxies. These non-gravitational processes are thought to raise the temperature of the IGM up to $`10^7`$K, as it seems required by the existence of an “entropy floor” in the scaling relation of groups of clusters (Ponman et al., 1999). In that case, several inter-galactic shocks would be suppressed and the $`\gamma `$-ray background flux generated there substantially reduced (Totani & Inoue, 2002). From our fig. 4 we infer a correction factor with respect to our estimate of order $`1/3`$. However, the actual role of such feed-back processes is still highly controversial (e.g., Kravtsov & Yepes, 2000, and references therein) and alternative scenarios in which the presence of cooling reduces dramatically the amount of required pre-heating are being investigated (Voigt & Bryan, 2001, and references therein).

According to our result in fig. 3, the the Gamma-ray Large Array Space Telescope should be able to measure the $`\gamma `$-ray flux due to IC emission from CR electrons accelerated at cluster accretion shocks (see also Miniati, 2002). In a previous work we explored with the same numerical technique adopted here the role of cluster accretion shocks as sites where the CR electrons that power radio relics (and not radio halos) are being produced (Miniati et al., 2001a). Our numerical model was able to account for several observed emission properties of radio relics (see also Enßlin et al., 1998; Roettiger et al., 1999, for direct observations related to this model). Interestingly the upper limit on the electron acceleration efficiency found here, corresponding to a parameter $`R_{e/p}2\times 10^2`$, is in line with the values adopted in that work in order to reproduce the measured radio emission for a number of sources. This number is also in agreement with recent observations of supernova remnants by Allen et al. (2001), who find $`R_{e/p}1/160`$. In any case, upcoming observations both in soft and hard $`\gamma `$-ray should allow us to put better constraints on such parameters for clusters accretion shocks.

Thus, with these constraints on the acceleration parameters the contribution to the CGB from cosmological CRs is estimated at the level of 20 % and most likely below 30 % of the measured flux. We point out that according to recent detailed calculations of Strong et al. (2000) the diffuse continuum $`\gamma `$-ray emission from the Galaxy could be significantly higher than previously estimated, thus lowering the actual level of the diffuse extra-galactic CGB. Likely, the next generation of high sensitivity $`\gamma `$-ray observatories will be resolve part of the CGB into a faint population of AGNs. In any case, the flux contributed by cosmological CRs will perhaps turn out the only truly diffuse extra-galactic $`\gamma `$-ray background component.

## 5 Conclusions

We have investigated numerically the contribution to the CGB from cosmological CRs. We carried out a simulation of structure formation in the canonical $`\mathrm{\Lambda }`$CDM universe and followed directly the acceleration, spatial transport and energy losses/gains of three different CR components: shock accelerated ions and electrons as well as secondary e± . Our results can be summarized as follows:

* Cosmological CRs produce a significant fraction of order of 20 % and no more than 30 % of the measure CGB. The computed flux is dominated by a leptonic component ($``$ 70 %) and accounts for a non negligible contribution from hadronic processes ($``$ 30 %).
* A higher flux, although in principle admissible, would be in disagreement with experimental upper limits set by EGRET on $`\gamma `$-ray emission from nearby clusters of galaxies.
* Based on the above considerations, we set an upper limit on the efficiency of conversion of shock ram pressure into energy of relativistic CR electrons of order $`1\times [\mathrm{log}(p_{max}/m_ec)/\mathrm{log}(4\times 10^7)]\%`$, where $`p_{max}`$ is the maximum energy of the accelerated electrons.
* IC $`\gamma `$-ray emission from clusters of galaxies could be of the same order of magnitude as $`\gamma `$-ray emission from $`\pi ^0`$-decay.

## Acknowledgments

This work was carried out at the Max Planck Institut für Astrophysik under the auspices of the European Commission for the ‘Physics of the Intergalactic Medium. I am indebted to O. Reimer for allowing me to use his unpublished results. I am grateful to A. Strong, O. Reimer, T. Enßlin and T. W. Jones for useful comments on the paper. The Max-Planck-Gesselschaft Rechenzentrum in Garching is acknowledged for precious support.


text:
# I Introduction

## I  Introduction

A few years ago, in Ref.  a method was devised to find alternative Lagrangians for the time dependent oscillator

$$\stackrel{..}{q}+\omega (t)^2q=0,$$
(1)
where $`q=q(t),`$ $`\omega (t)`$ is a given differentiable function, and the dot stands for time-derivative. The method, which is based on the concept of fouling transformation , is reviewed in Section II.

In this paper we study certain nonconventional quantum Hamiltonians corresponding to the classical fouled Hamiltonians associated with the Lagrangians derived in Ref. . This investigation is motivated by the fact that these nonconventional quantum Hamiltonians, having a polynomial structure in the operators $`a`$ and $`a^{},`$ may play an important role in the context of quantum optics especially in handling coherence and squeezing of multiphoton systems.

All the fouled Lagrangians of the hierarchy found in Ref.  lead to the same equation of motion (1) as it occurs for the conventional Lagrangian $`L_1=\frac{1}{2}(\stackrel{2}{\stackrel{.}{q}}\omega (t)^2q^2)`$ (see (5)). In Section III the fouled Hamiltonians $`K_\pm ,`$ related to the simplest fouled Lagrangians, $`L_2^{(1)}`$ and $`L_2^{(2)},`$ given by (23) and (24) are written down. By way of example, we have limited ourselves to consider the standard (harmonic) oscillator where $`\omega (t)\lambda `$ is a constant. Furthermore, it is shown that (at the classical level), as one expects, $`K_\pm `$ reproduce the same equation of motion coming from the conventional Hamiltonian $`H_0=\frac{1}{2}(p^2+\lambda ^2q^2`$ $`)`$. For brevity, we shall deal with $`K_+`$ only. It turns out that $`K_+`$ (see Section IV) can assume two possible forms, denoted by $`K_+^{(1)}`$ $`K_1`$ and $`K_+^{(2)}`$ $`K_2`$ (see (39) and (40)), which are independent but depend explicitly on time via the coefficients $`\mathrm{cos}\lambda t`$ and $`\mathrm{sin}\lambda t.`$ Consequently, $`K_1`$ and $`K_2`$ can be formally interpreted as the result of a rotation (of an angle $`\lambda t`$) of the quantities $`H_1=\sqrt{\lambda }(p^2+\lambda ^2q^2)q`$ and $`H_2=\frac{2}{3\sqrt{\lambda }}p^3,`$ which are not explicitly dependent on time. In other words, we have $`K_1^2+K_2^2=H_1^2+H_2^2,`$ and

$$\{K_1,K_2\}_{q,p}=\{H_1,H_2\}_{q,p}$$
(2)
(see Section IV), where the symbol $`\{,\}`$ denotes the Poisson bracket with respect to the canonical variables $`(q,p),`$ namely

$$\{A,B\}_{q,p}=\frac{A}{q}\frac{B}{p}\frac{A}{p}\frac{B}{q}.$$
(3)
In Section V the canonical quantization procedure based on a bosonic couple of annihilation and creation operators, $`a`$ and $`a^{}`$ (with $`[a,a^{}]=1),`$ is applied to the fouled Hamiltonians $`K_1`$ and $`K_2.`$ In such a way $`K_1`$ and $`K_2`$ are converted into a pair of Hamiltonian operators, say $`𝒦_1`$ and $`𝒦_2.`$ (Classically, as we know, $`K_1`$ and $`K_2`$ provide the same equation of motion (1)). Furthermore, since the canonical quantization procedure does affect only the cubic polynomial functions of $`(q,p),`$ namely $`H_1`$ and $`H_2,`$ here we limit ourselves to study the operators $`_1`$ and $`_2,`$ which are the quantized versions of $`H_1`$ and $`H_2,`$ respectively. The arising quantum models are investigated in Sections VI and VII.

Precisely, in Section VI we show that the quantum model coming from $`_1`$ gives rise to a linear second order eigenvalue equation of the Sturm-Liouville type (in the harmonic oscillator excitation number representation, or $`n`$rep).

An interesting feature of the operator $`_1`$ is that it is connected with an undetermined Hamburger moment problem . We show that this operator has deficiency indices (1,1) and allows a one-parameter family of self-adjoint extensions whose spectra are discrete and with no point in common. This goal has been achieved essentially by associating the operator $`_1`$ with a Jacobi matrix . A possible physical interpretation of the operator $`_1`$ is provided by fixing, among the infinite self-adjoint extensions, the extension corresponding to $`\epsilon =0,`$ where $`\epsilon `$ denotes the spectral parameter. In this case the eigenvalue equations for $`_1`$ can be solved exactly in all the following representations: the harmonic oscillator excitation number representation ($`n`$ rep), the coordinate space representation ($`q`$ rep), and the Fock-Bargmann holomorphic function representation ($`z`$ rep) . In all these representations, for $`\epsilon =0`$ square-integrable solutions of the eigenvalue equations are explicitly obtained.

In Section VII we deal with the eigenvalue equation for $`_2.`$ The analysis of this Hamiltonian is trivial. It is exactly solvable , its spectrum is on the whole line and the related eigenfunctions can be easily found. Nevertheless, $`_2`$ has been involved in the construction of two quantum models described by the operators $`_3`$ and $`_4`$ (see (57) and (58)), where $`_2`$ $`=`$ $`_3`$ $`+`$ $`_4.`$ One can attribute to $`_3`$ a physical meaning. Specifically, $`_3`$ can be interpreted as a special case of a class of Hamiltonians appearing in the higher order nonlinear optical processes . This Hamiltonian is involved in the construction of third power squeezed states .

On the other hand, as we can see in Section VII, the operator $`_4`$ is closely related to $`_1`$ so that the solutions of the corresponding eigenvalue equation can be derived from the solutions of the eigenvalue equation for $`_1`$ .

In Section VIII some constants of motion involving the operators $`_1`$,$`_2,_3,_4`$ and $`_5`$ $`a^3+a^3`$ are derived. Section IX contains a discussion on a possible algebraic framework which could be employed to analyze Eq. (69). This approach resorts to a quadratic algebra in terms of which the operator $`_1`$ can be naturally expressed. Finally, in Section X a few concluding remarks are presented, while in the Appendix a detail of a calculation is reported.

## II  <br>Fouled Lagrangians

We recall that a fouling transformation is a transformation under which the coordinates in configuration space are preserved:

$$Q=q,$$
(4)
while

$$P_n=P_n(q,p,t),$$
(5)
$`P_n`$ being a polynomial of degree $`n`$ in the variables $`q`$ and $`p=\frac{L_1}{\stackrel{.}{q}}=\stackrel{.}{q},`$ i.e.

$$P_n=\underset{0}{\overset{n}{}}a_jp^{nj}q^j,$$
(6)
where

$$L_1=\frac{1}{2}(\stackrel{.}{q^2}\omega (t)^2q^2)$$
(7)
is the conventional Lagrangian, and $`a_j=a_j(t)`$ are (real) time-dependent coefficients.

In  it was proven that the function $`L_n`$ $`=L_n(q,\stackrel{.}{q},t)`$ expressed by

$$L_n=\underset{j=0}{\overset{n}{}}\frac{1}{nj+1}a_j\stackrel{nj+1}{\stackrel{.}{q}}q^j+(\underset{n}{\overset{.}{a}}\omega ^2a_{n1})\frac{q^{n+1}}{n+1}$$
(8)
satisfies the equations

$$\frac{L_n}{\stackrel{.}{q}}=P_n,\frac{L_n}{q}=\underset{n}{\overset{.}{P}}.$$
(9)
The compatibility condition for these equations provides

$$\underset{o}{\overset{.}{a}}=\frac{n1}{n}a_1,\underset{j}{\overset{.}{a}}=(nj+1)\omega ^2a_{j1}(j+1)\frac{nj1}{nj}a_{j+1},$$
(10)
with $`j=1,2,\mathrm{},n1.`$ Furthermore, we have

$$\frac{d}{dt}\frac{L_n}{\stackrel{.}{q}}\frac{L_n}{q}=I_n(t)(\stackrel{..}{q}+\omega (t)^2q),$$
(11)
where $`I_n(t)`$ is a time-dependent constant of motion, viz. $`\frac{d}{dt}I_n(t)`$ $`=0,`$ given by

$$I_n(t)=\underset{j=0}{\overset{n}{}}(nj)a_jp^{nj1}q^j=\frac{^2L_n}{\stackrel{2}{\stackrel{.}{q}}},$$
(12)
for any $`n1.`$

From Eq. (9) we deduce that corresponding to the solution $`q`$ of the generalized oscillator (1), $`L_n`$ satisfies the Euler-Lagrange equation. We notice that for $`n=1`$, Eq. (6) gives the conventional Lagrangian (5) (with $`a_0=1`$ and $`a_1=0).`$ The related invariant is $`I_1=\frac{^2L_1}{\stackrel{2}{\stackrel{.}{q}}}=1.`$

Hereafter, by way of example, we are interested in the case $`n=2.`$ Then, from (6) and (10) we obtain

$$L_2=\frac{1}{3}a_0\stackrel{3}{\stackrel{.}{q}}+\frac{1}{2}a_1\stackrel{2}{\stackrel{.}{q}}q+a_2\stackrel{.}{q}q^2+\frac{1}{3}(\stackrel{.}{a_2}\omega ^2a_1)q^3,$$
(13)
$$I_2=2a_0p+a_1q,$$
(14)
where (see (8))

$$\underset{0}{\overset{.}{a}}=\frac{1}{2}a_1,\underset{1}{\overset{.}{a}}=2\omega ^2(t)a_0,$$
(15)
namely

$$\underset{0}{\overset{..}{a}}+\omega ^2(t)a_0=0.$$
(16)
The general solution of Eq. (14) can be written in the form

$$a_0=\sqrt{2}\sigma (c_1\mathrm{cos}\frac{\theta }{2}+c_2\mathrm{sin}\frac{\theta }{2}),$$
(17)
where $`c_1,c_2`$ are constants, and $`\sigma ,\theta `$ are defined by

$$\stackrel{..}{\sigma }+\omega ^2(t)\sigma =\frac{1}{4\sigma ^3},$$
(18)
and

$$\stackrel{.}{\theta }=\frac{1}{\sigma ^2}.$$
(19)
In the following we shall limit ourselves to the choice $`\omega (t)=\lambda =const.`$ Consequently, Eqs. (16) and (17) admit the solution

$$\sigma =\frac{1}{\sqrt{2\lambda }}$$
(20)
and

$$\theta =2\lambda t.$$
(21)
Equation (15) has the two independent solutions

$$a_0^{(1)}=\frac{1}{\sqrt{\lambda }}\mathrm{cos}\lambda t,a_0^{(2)}=\frac{1}{\sqrt{\lambda }}\mathrm{sin}\lambda t,$$
(22)
while the corresponding values $`a_1^{(1)},`$ $`a_1^{(2)}`$ are (see (13))

$$a_1^{(1)}=2\sqrt{\lambda }\mathrm{sin}\lambda t,a_1^{(2)}=2\sqrt{\lambda }\mathrm{cos}\lambda t.$$
(23)
On the other hand, the expressions of $`a_2^{(1)},`$ $`a_2^{(2)}`$ turn out to be (see )

$$a_2^{(1)}=\lambda ^{\frac{3}{2}}\mathrm{cos}\lambda t,a_2^{(2)}=\lambda ^{\frac{3}{2}}\mathrm{sin}\lambda t.$$
(24)
Now, using (20), (21) and (22), from (11) we obtain the two alternative fouled Lagrangians

$$L_2^{(1)}=\frac{1}{3\sqrt{\lambda }}\mathrm{cos}\lambda t\stackrel{3}{\stackrel{.}{q}}+\sqrt{\lambda }\stackrel{}{\mathrm{sin}\lambda t\stackrel{.}{q^2}}q+\lambda ^{^{\frac{3}{2}}}\mathrm{cos}\lambda t\stackrel{.}{q}q^2\lambda ^{\frac{5}{2}}\mathrm{sin}\lambda tq^3,$$
(25)
$$L_2^{(2)}=\frac{1}{3\sqrt{\lambda }}\mathrm{sin}\lambda t\stackrel{3}{\stackrel{.}{q}}\sqrt{\lambda }\stackrel{}{\mathrm{cos}\lambda t\stackrel{.}{q^2}}q+\lambda ^{\frac{3}{2}}\mathrm{sin}\lambda t\stackrel{.}{q}q^2+\lambda ^{\frac{5}{2}}\mathrm{cos}\lambda tq^3,$$
(26)
which furnish the two independent invariants

$$I_2^{(1)}=\frac{^2L_2^{(1)}}{\stackrel{2}{\stackrel{.}{q}}}=\frac{2}{\sqrt{\lambda }}\mathrm{cos}\lambda t\stackrel{.}{q}+2\sqrt{\lambda }\mathrm{sin}\lambda tq,$$
(27)
$$I_2^{(2)}=\frac{^2L_2^{(2)}}{\stackrel{2}{\stackrel{.}{q}}}=\frac{2}{\sqrt{\lambda }}\mathrm{sin}\lambda t\stackrel{.}{q}2\sqrt{\lambda }\mathrm{cos}\lambda tq,$$
(28)
respectively (see (10)).

Equations (25) and (26) provide the general solution $`q`$ and the momentum $`p=\stackrel{.}{q}`$ of the harmonic oscillator ($`\omega \lambda ),`$ i.e.

$$q=\frac{1}{2\sqrt{\lambda }}(I_2^{(1)}\mathrm{sin}\lambda tI_2^{(2)}\mathrm{cos}\lambda t),$$
(29)
$$p=\stackrel{.}{q}=\frac{\sqrt{\lambda }}{2}(I_2^{(1)}\mathrm{cos}\lambda t+I_2^{(2)}\mathrm{sin}\lambda t),$$
(30)
from which, as one expects,

$$H_0=\frac{1}{2}(p^2+\lambda ^2q^2)=\frac{\lambda }{8}[(I_2^{(1)})^2+(I_2^{(2)})^2]=const,$$
(31)
where $`H_0`$ is the conventional Hamiltonian.

## III Fouled Hamiltonians

At this point, we build up the Hamiltonian corresponding to the Lagrangian $`L_2`$ given by (11). Following the procedure of Ref. , we get

$$K_\pm =\frac{a_1}{2a_0}qP\pm \frac{2}{3}a_0[(\frac{a_1^2}{4a_0^2}\frac{a_2}{a_0})q^2+\frac{P}{a_0}]^{\frac{3}{2}}+[\frac{a_1}{6a_0}(3a_2\frac{a_1^2}{2a_0})\frac{1}{3}(\stackrel{.}{a_2}\lambda ^2a_1)]q^3,$$
(32)
where

$$P=a_0p^2+a_1qp+a_2q^2.$$
(33)
For the sake of definiteness, later we shall study the Hamiltonian $`K_+`$ only.

The Hamilton equations for $`K_+`$ read

$$\stackrel{.}{q}=\frac{a_1}{2a_0}q+\frac{I_2}{2a_0},$$
(34)
$$\stackrel{.}{P}=\frac{a_1}{2a_0}P\frac{a_1^24\lambda ^2a_0^2}{4a_0^2}I_2q3a_1(\lambda ^2\frac{a_1^2}{12a_0^2})q^2,$$
(35)
where the relation

$$(\frac{a_1^2}{4a_0^2}\frac{a_2}{a_0})q^2+\frac{P}{a_0}=(\frac{I_2}{2a_0})^2$$
(36)
with $`a_2=\lambda ^2a_0`$ has been used (see ).

By expliciting Eqs. (31) and (32) with the help of (12) and (30), we arrive at the equation

$$I_2(\stackrel{..}{q}+\lambda ^2q)=0,$$
(37)
which tells us that the fouled Hamiltonian $`K_+,`$ as it occurs for the related fouled Lagrangian $`L_2`$ (see (9)), gives rise to the same equation of motion emerging from the conventional Hamiltonian (29).

By virtue of (33), the Hamiltonians (30) take the form

$$K_+=\frac{1}{2}a_1(p^2+\lambda ^2q^2)q+\frac{2}{3}a_0p^3,$$
(38)
$$K_{}=\frac{3}{2}a_1p^2q\frac{a_1^2}{a_0}q^2p+(\frac{a_1}{2}\lambda ^2\frac{a_1^3}{6a_0^2})q^3\frac{2}{3}a_0p^3.$$
(39)
One can see that

$$K_\pm =PpL_2,$$
(40)
where $`P`$ and $`L_2`$ are given by (31) and (11).

## IV The fouled Hamiltonians $`K_+^{(1)}`$ and $`K_{}^{(2)}`$

By using (20), (21) and (22), from (36) we derive the pair of (explicitly time-dependent) Hamiltonians

$$K_+^{(1)}K_1=H_1\mathrm{sin}\lambda t+H_2\mathrm{cos}\lambda t,$$
(41)
$$K_+^{(2)}K_2=H_1\mathrm{cos}\lambda t+H_2\mathrm{sin}\lambda t,$$
(42)
where $`H_1`$ and $`H_2`$ are defined by

$$H_1=\sqrt{\lambda }(p^2+\lambda ^2q^2)q=2\sqrt{\lambda }H_0q,$$
(43)
$$H_2=\frac{2}{3\sqrt{\lambda }}p^3.$$
(44)
One can check straightforwardly that the following evolution equations

$$\stackrel{.}{K_1}=\{K_1,H_0\}_{q,p}+\frac{K_1}{t},$$
(45)
$$\stackrel{.}{K_2}=\{K_2,H_0\}_{q,p}+\frac{K_2}{t},$$
(46)
$$\stackrel{.}{H_1}=\{H_1,H_0\}_{q,p},$$
(47)
$$\stackrel{.}{H_2}=\{H_2,H_0\}_{q,p},$$
(48)
hold, where $`H_0`$ is given by (29) and

$$\{H_1,H_0\}_{q,p}=2\sqrt{\lambda }pH_0,$$
(49)
$$\{H_2,H_0\}_{q,p}=2\frac{\lambda ^2}{\sqrt{\lambda }}p^2q.$$
(50)
We also have that the Poisson brackets between $`K_1,K_2`$ and $`H_1,H_2`$ concide, i.e.

$$\{K_1,K_2\}_{q,p}=\{H_1,H_2\}_{q,p}.$$
(51)
This is a direct consequence of the rotation form of the transformations (39) and (40).

## V Quantization

Hereafter, we shall put for simplicity $`\lambda =\mathrm{}=1.`$

In order to quantize the fouled Hamiltonians (39) and (40), let us introduce the operators

$$\stackrel{}{q}=\frac{1}{\sqrt{2}}(a+a^{}),\stackrel{}{p}=\frac{i}{\sqrt{2}}(aa^{}),$$
(52)
where $`a`$ and $`a^{}`$ denote a (boson) annihilation and a creation operator, respectively.

By means of (50), we can write the operators $`_1`$and $`_2`$ corresponding to the classical functions (41) and (42). We obtain

$$_1=\sqrt{2}(a^2a+a^{}a^2+a^{}+a),$$
(53)
$$_2=\frac{i}{\sqrt{2}}[\frac{1}{3}(a^3a^3)+(a^2aa^{}a^2+a^{}a)],$$
(54)
with the help of the commutation rule

$$[a,a^{}]=1$$
(55)
or, in terms of the Heisenberg commutation relation: $`[\stackrel{}{q},\stackrel{}{p}]=i.`$

Now let us use the representation $`\stackrel{}{q}=x,`$ $`\stackrel{}{p}=i\frac{d}{dx},`$ so that the operators $`a`$ and $`a^{}`$ can be written as (see (50))

$$a=\frac{1}{\sqrt{2}}(x+\frac{d}{dx}),a^{}=\frac{1}{\sqrt{2}}(x\frac{d}{dx}).$$
(56)
In terms of these quantities, $`_1`$and $`_2`$ take the forms

$$_1=(x\frac{d^2}{dx^2}+\frac{d}{dx})+x^3,$$
(57)
and

$$_2=\frac{2i}{3}\frac{d^3}{dx^3}.$$
(58)
For later convenience, we shall report also in the representation (54) the following operators appearing in (52):

$$_3\frac{i}{\sqrt{2}}\frac{1}{3}(a^3a^3)=\frac{i}{2}[\frac{1}{3}\frac{d^3}{dx^3}+(x^2\frac{d}{dx}+x)],$$
(59)
$$_4\frac{i}{\sqrt{2}}(a^2aa^{}a^2+a^{}a)=\frac{i}{2}[\frac{d^3}{dx^3}(x^2\frac{d}{dx}+x)].$$
(60)
In the next Section, we shall study the operators $`_1`$and $`_2`$ by dealing with the corresponding eigenvalue problems.

## VI  Self-adjoint extensions of the operator $`_1`$

In order to clarify the quantum-mechanical meaning of the Hamiltonian $`_1,`$ it is crucial to establish whether $`_1`$ enjoys self-adjoint type properties. In doing so, first let us recall that one can provide different representations for the operator (51), which correspond to various forms of the related eigenvalue equations. We shall consider the following representations: the harmonic oscillator excitation number representation ($`n`$ rep), the coordinate representation ($`q`$ rep), and the Fock-Bargmann holomorphic function reprentation ($`z`$ rep).

For reader’s convenience, we shall summarize below the main properties of these representations .

Let us $`𝐇`$ denote the Hilbert space where the operators $`\stackrel{}{q},`$ $`\stackrel{}{p}`$ and $`a^{},`$ $`a`$ act. In the $`n`$ rep, the vectors $`\{n>\}`$ form a basis in $`𝐇`$. The following relations

$$a=\sqrt{n}n1>,a^{}n>=\sqrt{n+1}n+1>,a^{}an>=nn>$$
(61)
hold.

On the other hand, in the $`q`$ rep a vector $`\psi >`$ belonging to the space $`𝐇`$ is represented by a coordinate function $`<q\psi >=\psi (q)`$ which is square-integrable:

$$_{\mathrm{}}^+\mathrm{}\psi (q)^2𝑑q<\mathrm{}.$$
(62)
The basic vector $`n>`$ is described by the function

$$<qn>=\phi _n(q)=N_nH_n(q)\mathrm{exp}(\frac{q^2}{2}),$$
(63)
where $`N_n=(\sqrt{\pi }2^nn!)^{\frac{1}{2}}`$ and $`H_n(q)`$ is the Hermite polynomial of degree $`n.`$ In the $`q`$ rep, formulas (59) become the standard recursion relations for the Hermite polynomials.

In order to introduce the Fock-Bargmann representation (or $`z`$ rep), let $`\psi >`$ be an arbitrary normalized vector in $`𝐇,`$ namely

$$\psi >=\underset{n=0}{\overset{\mathrm{}}{}}c_nn>,$$
(64)
with $`<\psi \psi >=_{n=0}^{\mathrm{}}c_n^2=1.`$ Furthermore, taking into account the Glauber form

$$z>=\mathrm{exp}(\frac{z^2}{2})\underset{n=0}{\overset{\mathrm{}}{}}\frac{z^n}{\sqrt{n!}}n>,$$
(65)
the state $`\psi >`$ is completely determined by

$$<z\psi >=\mathrm{exp}(\frac{z^2}{2})\psi (\stackrel{}{z}),$$
(66)
where

$$\psi (z)=\underset{n=0}{\overset{\mathrm{}}{}}c_nu_n(z),$$
(67)
with $`u_n(z)=\frac{z^n}{\sqrt{n!}}.`$ Owing to the condition $`_{n=0}^{\mathrm{}}c_n^2=1,`$ the series in (65) converges uniformly in any compact domain of the complex $`z`$ plane. Consequently, $`\psi (z)`$ turns out to be an entire holomorphic function in the $`z`$ plane, and

$$\psi ^2=\mathrm{exp}(z^2)\psi (z)^2𝑑\mu (z)<\mathrm{},$$
(68)
where $`d\mu (z)=\pi ^1dxdy,`$ $`z=x+iy`$ .

The scalar product of two entire functions $`\psi _1(z)`$ and $`\psi _2(z)`$ obeying the condition (66) is given by

$$<\psi _1\psi _2>=\mathrm{exp}(z^2)\underset{1}{\overset{}{\psi }}(z)\psi _2(z)d\mu (z).$$
(69)
As it was proven by Bargmann , the Fock-Bargmann representation space with a scalar product provided by (67), is really a Hilbert space. We observe also that in the $`z`$rep, the operator solution for the commutation relation $`[a,a^{}]=1`$  is

$$a\frac{d}{dz},a^{}z.$$
(70)
Now let us consider the eigenvalue equation

$$(a^2a+a^{}a^2+a+a^{})\chi >=\epsilon \chi >.$$
(71)
Starting from the $`n`$ rep and following the lines of Ref. , let us put

$$\chi >=\underset{n=0}{}f_n(\epsilon )n>$$
(72)
into Eq. (69). By using Eqs. (59), after simple calculations we obtain the recursion formula

$$(n+1)^{\frac{3}{2}}f_{n+1}(\epsilon )\epsilon f_n(\epsilon )+n^{\frac{3}{2}}f_{n1}(\epsilon )=0$$
(73)
for $`n1,`$ and

$$f_1(\epsilon )=\epsilon f_0(\epsilon )$$
(74)
(the boundary condition).

We remark that the sequence $`\{f_n\}`$ is such that the series $`_{n=0}f_n^2`$ converges, i.e. the sequence $`\{f_n\}(f_0,f_1,f_2,\mathrm{})`$ belongs to the Hilbert space $`l^2.`$

To show this, first let us consider the case $`\epsilon =0.`$ Then, Eqs. (71)-(72) provide

$$f_{2n}=(1)^n[\frac{(2n1)!!}{(2n)!!}]^{\frac{3}{2}}f_0,$$
(75)
the odd terms being zero. Thus, from the asymptotic formula of Gamma function (, p. 257). we deduce that for $`n\mathrm{},`$ $`f_{2n}^2`$ behaves as $`n^{\frac{3}{2}},`$ so that $`_{n=0}f_{2n}^2<\mathrm{}.`$ In general, i.e. for $`\epsilon 0,`$ Eq. (71) tells us that for great values of $`n,`$ the sum of the first and the last term is leading with respect to the second term. This allows us to see easily that both even and odd terms, $`f_{2n}^2`$ and $`f_{2n+1}^2,`$ behave asymptotically as $`n^{\frac{3}{2}}.`$ To conclude, the series $`_{n=0}f_n^2`$ is convergent, namely the sequence $`\{f_n\}`$ belongs to the Hilbert space $`l^2`$ for any value of the spectral parameter $`\epsilon .`$ Then, the equation $`_1^{}\chi >=\epsilon \chi >,`$ for $`\mathrm{}\lambda 0,`$ has nontrivial solutions (, p. 140).

Now, by introducing the notation $`b_n=(n+1)^{\frac{3}{2}},`$ we see that it is possible to associate with the difference equation (71) the (infinite) Jacobi matrix

$$A=\left(\begin{array}{ccccc}0& b_0& 0& 0& \mathrm{}\\ b_0& 0& b_1& 0& \mathrm{}\\ 0& b_1& 0& b_2& \mathrm{}\\ 0& 0& b_2& 0& \mathrm{}\\ ..& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right)$$
(76)
so that Eq. (69) is equivalent to the eigenvalue equation $`Af=\epsilon f,`$ with $`f=(f_0,f_1,`$ $`\mathrm{})^T.`$

The Jacobi matrix (74) plays a crucial role in the study of the Hamburger moment problem (see, for example, ). Precisely, let us consider the moments

$$s_n=_{\mathrm{}}^+\mathrm{}x^n𝑑\sigma (x),n=0,1,2,\mathrm{},$$
(77)
where $`\sigma `$ denotes a (positive) measure on $``$ (, p. 145). The Hamburger moment problem is to determine conditions on a sequence of real numbers $`\stackrel{}{\{s_n\}_{n=0}^{\mathrm{}}}`$ , so that there exists a measure satisfying (75). One can show that a sequence of real numbers $`\{s_n\}`$ are the moments of a positive measure on $``$ if and only if for all $`N`$ and all $`\alpha _0,\alpha _1,\mathrm{},\alpha _N`$ $`𝒞`$, one has

$$\underset{n,m=0}{\overset{N}{}}\underset{n}{\overset{}{\alpha }}\alpha _ms_{n+m}0.$$
(78)
From the Jacobi matrix (74) we get the limitations

$$\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{b_n}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{1}{(n+1)^{\frac{3}{2}}}<\mathrm{},b_{n1}b_{n+1}<b_n^2.$$
(79)
Consequently, the Jacobi matrix (74) belongs to the $`type`$ $`C`$ (limit circle case), and corresponds to an undetermined Hamburger moment problem (). So, the properties of the operator $`_1`$ $``$ $`a^2a+a^{}a^2+a+a^{}`$ are similar to the properties of the operator $`a^k+a^k`$ ($`k=3)`$ discussed by Nagel . In other words, the operator $`_1`$ has deficiency indices (1,1) and allows a one-parameter family of self-adjoint extensions, each having a purely discrete spectrum on the real line . The spectra of two different extensions turn out to have no point in common ( , p. 152). We have that different self-adjoint extensions correspond to different dynamics .

Since every self-adjoint extension of $`_1`$ has a discrete spectrum on the real line, and taking one eigenvalue determines the corresponding extension uniquely, let us choose $`\epsilon =0.`$

In this case the eigenvalue equation underlying the operator $`_1`$ can be solved exactly in all the representations mentioned at the beginning of this Section.

To show this, let us deal with the $`q`$ rep. With the help of (54), the eigenvalue equation (69) reads

$$x\varphi ^{\prime \prime }+\varphi ^{}+(\sqrt{2}\epsilon x^3)\varphi =0,$$
(80)
where $`\varphi ^{}\frac{d}{dx}\varphi `$.

Equation (78) can be written as the Sturm-Liouville equation (, p. 59)

$$[\varphi (x)]=\sqrt{2}\epsilon \varphi (x),$$
(81)
where $``$ denotes the Sturm-Liouville operator

$$=\frac{d}{dx}[x\frac{d}{dx}]x^3.$$
(82)
By using the transformation

$$\varphi (x)=\mathrm{exp}(\frac{x^2}{2})\psi (x),$$
(83)
Equation (78) becomes

$$x\psi ^{\prime \prime }+(12x^2)\psi ^{}+(\sqrt{2}\epsilon 2x)\psi =0.$$
(84)
This equation is satisfied by

$$\psi (x)=\underset{n=0}{\overset{\mathrm{}}{}}f_nN_nH_n(x),$$
(85)
where the coefficients $`f_n`$ fulfil the recursion relations (71) and (72). We have already shown that $`\{f_n\}`$ $`l^2.`$ Then, the function $`\psi (x)`$ belongs to the Hilbert space $`L_{e^{x^2}}^2(\mathrm{},\mathrm{})`$ ($`\mathrm{exp}(x^2)`$ is the weight function). We point out that generally the relation (83) is not valid in the pointwise sense, but it holds in accordance with the metric of $`L_{e^{x^2}}^2(\mathrm{},\mathrm{}),`$ namely

$$\underset{n\mathrm{}}{lim}_{\mathrm{}}^+\mathrm{}\psi (x)\underset{k=0}{\overset{n}{}}f_kN_kH_k(x)^2\mathrm{exp}(x^2)𝑑x=0.$$
(86)
For $`\epsilon =0,`$ via the change of variable $`\xi =x^2`$ Eq. (82) can be written as a special case of a Kummer equation, whose independent solutions are $`M(\frac{1}{2},1;x^2)`$ and $`U(\frac{1}{2},1;x^2)`$ ( , p. 504). In the case the solution $`\psi (x)`$ of Eq. (82) with the property $`\psi (x)`$ $`L_{e^{x^2}}^2(\mathrm{},\mathrm{})`$ is given by

$$\psi (x)=cU(\frac{1}{2},1;x^2)=\underset{n=0}{\overset{\mathrm{}}{}}f_{2n}N_{2n}H_{2n}(x),$$
(87)
where the constant $`c`$ is such that $`c\pi ^{\frac{3}{4}}=f_0,`$ and

$$f_{2n}=f_0\pi ^{\frac{3}{4}}N_{2n}_{\mathrm{}}^+\mathrm{}U(\frac{1}{2},1;x^2)H_{2n}(x)\mathrm{exp}(x^2)𝑑x.$$
(88)
One can easily check that (86) is satisfied for any $`n𝒩`$.

Finally, in the $`z`$ rep, i.e. for $`\chi >`$ $`\chi ,`$ $`a^{}z,`$ $`a\frac{d}{dz},`$ Eq. (69) gives

$$z\chi _{zz}+(1+z^2)\chi _z+(z\epsilon )\chi =0,$$
(89)
where

$$\chi =\underset{n=0}{\overset{\mathrm{}}{}}f_n\frac{z^n}{\sqrt{n!}}$$
(90)
and $`f_n`$ satisfies the recursion relations (71) and (72). For $`\epsilon =0,`$ the eigenvalue equation (87) as well can be exactly solved. In fact, by setting $`z^2=y,`$ $`\zeta =\frac{y}{2},`$ this equation becomes a special case of the Kummer equation.

Therefore, in the $`z`$ rep, where the eigenfunction should be a holomorphic (and normalizable) function in the whole $`z`$ plane, one has the solution $`M(\frac{1}{2},1;\frac{z^2}{2}).`$ The other solution, namely the Kummer function $`U(\frac{1}{2},1;\frac{z^2}{2})`$, is not holomorphic at $`z=0.`$

To conclude this Section, we observe that the solutions of the eigenvalue equation for $`_1`$ in the case $`\epsilon =0`$ can be also found from Eqs. (71) and (72), obtained within the $`n`$ rep, by means of the standard integral representations of the confluent hypergeometric functions. An example of this procedure is displayed in .

## VII The $`_2`$, $`_3`$, $`_4`$ models

The eigenvalue problem for the operator $`_2`$ can be written as

$$\psi =\epsilon \psi ,$$
(91)
where $`=\frac{2i}{3}D_x^3`$ $`(D_x=\frac{d}{dx}).`$ The study of Eq. (89) is trivial. It is exactly solvable , its spectrum is on the whole line, and the related eigenfunctions are of the exponential type.

The operator $`_3`$ (see (57)) belongs to the class of Hamiltonians

$$=i\kappa _n(a^na^n)$$
(92)
appearing in the higher order nonlinear optical processes. In particular, (57) describes a subharmonic generation process, in which a photon from a strong pump beam produces $`n`$ photons of the signal beam in a nonlinear medium . The constant $`\kappa _n`$ is related to the $`n`$th nonlinear susceptibility coefficient and to the amplitude of the pump field, while $`a`$ and $`a^{}`$ are the annihilation and the creation operators for the signal field. In this context, the evolution of an arbitrary initial state $`\mathrm{\Psi }(0)>`$ of the signal field to the state $`\mathrm{\Psi }(t)>`$ is governed by

$$\mathrm{\Psi }(t)>=\mathrm{exp}[\kappa _nt(a^na^n)]\mathrm{\Psi }(0)>.$$
(93)
The squeezing of this state was examined by Hillery, Zubairy and Wódkiewicz . They showed that to any order in the coupling constant $`\kappa _n,`$ the vacuum state is not squeezed in the higher order nonlinear optical processes $`(n3).`$

This important result stimulated the analysis of $`n`$th power squeezed states. Interesting (and, generally, not yet completely explored) questions arise in connection with this argument. Some of them are discussed in  and references therein.

Now let us make some comments about the Hamiltonian $`_4.`$ This operator is closely related to $`_1.`$ This can be seen by means of the phase transformation $`a^{}=ia,`$ $`a^{}=ia^{},`$ so that $`_4`$ takes the form

$$_4=\frac{1}{\sqrt{2}}(a^2a^{}+a^{}a^2+a^{}+a^{}).$$
(94)
In other words, one has $`_4=\frac{1}{2}_1`$ (in terms of the primed operators). This corresponds to pick up $`q^{}=p`$ and $`p^{}=q.`$ In such a way $`_4`$ turns out to be the inverse Fourier transform of (55). Therefore, the solutions of the eigenvalue equation for the Hamiltonian operator $`_4`$ can be derived from the solutions of the eigenvalue equation for the Hamiltonian operator $`_1`$ (see (69)).

The eigenvalue equation for $`_1`$ ( $`_4`$ ) can be investigated by means of the algebraic approach outlined in Section IX.

## VIII  Equations and constants of motion related to the operators $`_j`$

The equations of motion for the Hamiltonians $`_j`$ $`(j=1,2,3,4)`$ arise immediately by using the Heisenberg representation. In other words, by putting $`a(t)=a(0)\mathrm{exp}(it)`$ and $`a^{}(t)=a^{}(0)\mathrm{exp}(it)`$ in the expressions (51) and (58), we easily find (as one expects) that $`_1`$ and $`_4`$ sastisfy the same equation of motion (i.e., the equation for the harmonic oscillator of frequency $`\lambda =1`$):

$$\frac{d^2}{dt^2}_1+_1=0,\frac{d^2}{dt^2}_4+_4=0.$$
(95)
On the other hand, for the operators $`_3=\frac{i}{\sqrt{2}}\frac{1}{3}(a^3a^3)`$ (see (57)) and $`_5=\frac{1}{\sqrt{2}}\frac{1}{3}(a^3+a^3),`$ the same considerations made for $`_1`$ and $`_4`$ in Section VII hold. One has that $`_3`$ and $`_5`$ obey the same equation of motion (i.e., the equation for the harmonic oscillator of frequency $`3`$):

$$\frac{d^2}{dt^2}_3+9_3=0,\frac{d^2}{dt^2}_5+9_5=0.$$
(96)
Some comments on Eqs. (93) and (94) are presented in Section X.

At this point, we observe that in addition to the constants of motion

$$\stackrel{}{q}=e^{it\stackrel{}{H_0}}\stackrel{}{q}e^{it\stackrel{}{H_0}}=\stackrel{}{q}\mathrm{cos}t\stackrel{}{p}\mathrm{sin}t,$$
(97)
$$\stackrel{}{p}=e^{it\stackrel{}{H_0}}\stackrel{}{p}e^{it\stackrel{}{H_0}}=\stackrel{}{q}\mathrm{sin}t+\stackrel{}{p}\mathrm{cos}t,$$
(98)
where $`[\stackrel{}{q},\stackrel{}{p}]=i,`$ $`\stackrel{}{H_0}=\frac{1}{2}(\stackrel{2}{\stackrel{}{p}}+\stackrel{2}{\stackrel{}{q}})`$ and $`\stackrel{}{q}=\frac{1}{2}I_2^{(2)},`$ $`\stackrel{}{p}=\frac{1}{2}I_2^{(1)},`$ we obtain also the following set of constants:

$$\stackrel{}{H_1}=_1\mathrm{cos}t2_4\mathrm{sin}t,$$
(99)
$$\stackrel{}{H_4}=\frac{1}{2}_1\mathrm{sin}t+_4\mathrm{cos}t,$$
(100)
$$\underset{3}{\overset{}{H}}=_3\mathrm{cos}3t_5\mathrm{sin}3t,$$
(101)
$$\underset{5}{\overset{}{H}}=_3\mathrm{sin}3t+_5\mathrm{cos}3t.$$
(102)
In terms of the arbitrary constants $`\underset{j}{\overset{}{H}}(j=1,4,3,5),`$ one can express the general solutions of Eqs. (93) and (94), which can be easily found by inverting the transformations (97), (98) and (99), (100), respectively.

## IX A possible algebraic framework for the study of Eq. (69)

For $`\epsilon 0,`$ Eq. (69) could be investigated following different analytical techniques. Among these, an important role is played by the algebraic approach, which allows one to express Eq. (69) in terms of the generators of a quadratic algebra .

As an example of a quadratic algebra, we can consider a nonlinear deformation of the (linear) $`su(1,1)`$ algebra, defined by the relations

$$[J_0,J_\pm ]=\pm J_\pm ,$$
(103)
$$[J_+,J_{}]=P(J_0),$$
(104)
where the Jacobi identity holds. $`J_\pm `$ are the ladder operators, and $`P(J_0)`$ is a second-degree polynomial function of the diagonal operator $`J_0.`$ $`P(J_0)`$ can be written in the form

$$P(J_0)=\alpha _1+\alpha _2J_0+\alpha _3J_0^2,$$
(105)
where $`\alpha _i`$ ($`i=1,2,3)`$ are arbitrary coefficients ($`\alpha _30).`$

Now let us introduce the (bosonic) realization

$$J_0=a^{}a+\frac{1}{2},$$
(106)
$$J_{}=(a^{}a^2+a),$$
(107)
$$J_+=a^2a+a^{},$$
(108)
with $`J_+=(J_{})^{}.`$ Then, the operators (104)–(106) turn out to obey the commutation relations

$$[J_0,J_\pm ]=\pm J_\pm ,[J_+,J_{}]=\frac{1}{4}+3J_0^2.$$
(109)
We remark that the quadratic algebra (107) is a finite $`W_3^{(2)}`$algebra , which corresponds to choose $`\alpha _1=\frac{1}{4},`$ $`\alpha _2=0,`$ and $`\alpha _3=3`$ in (103).

By virtue of (104)–(106), Eq. (69) can be re-expressed as

$$\mathrm{\Omega }(\epsilon )\chi >=0,$$
(110)
with

$$\mathrm{\Omega }(\epsilon )=J_+J_{}\epsilon _.$$
(111)
Therefore, the eigenvalue problem (69) can be formulated in terms of the generators $`J_\pm `$ of a quadratic algebra of the type $`W_3^{(2)}.`$ The Casimir operator of the quadratic algebra (107) is considered in the Appendix.

## X  Conclusions

In the context of the existence of an infinite set of fouled Lagrangians and Hamiltonians for the generalized (time-dependent) oscillator, we have studied a special case where the frequency of the oscillator is assumed to be a constant $`\lambda `$ (harmonic oscillator). By way of example, we have considered a pair of independent fouled Lagrangians $`(L_2^{(1)}`$, $`L_2^{(2)})`$ (see (23)-(24)) and, in correspondence, a pair of fouled independent Hamiltonians $`(K_1,`$ $`K_2)`$ (see (39)-(40)) which lead, at the classical level, to the same equation of motion provided by the conventional Lagrangian and Hamiltonian. Both $`(L_2^{(1)}`$, $`L_2^{(2)})`$ and $`(K_1,`$ $`K_2)`$ are explicitly time-dependent.The method followed to find these alternative Lagrangians and Hamiltonians implies the construction of two independent invariants (constants of motion), $`I_2^{(1)}`$ and $`I_2^{(2)};`$ in terms of them the canonical variables $`q`$ and $`p`$ can be expressed. These invariants, which are of the Nöther type , are connected with the ”quadrature-phase amplitudes” appearing in the problem of generation of squeezed states in certain optical devices .

In this paper we have focused our attention mainly on the quantized version of the fouled Hamiltonians $`(K_1,`$ $`K_2).`$ Since $`K_1`$ and $`K_2`$ are related to the cubic polynomials $`H_1=\sqrt{\lambda }(p^2`$ $`+\lambda ^2q)q`$ and $`H_2=\frac{2}{3\sqrt{\lambda }}p^3`$ by a formal rotation, the canonical quantization prescription affects essentially $`H_1`$ and $`H_2.`$ Our purpose has been to study, at the quantum level, the Hamiltonians $`_1`$ and $`_2`$ (see (51), (52)) corresponding to $`H_1`$ and $`H_2`$ . The model described by $`_2`$ is associated with an eigenvalue problem given by a (linear) third order differential equation, with constant coefficients, which is exactly solvable . Furthermore, we have put $`_2=_3+_4`$ , where $`_4`$ (see (58)) turns out to be closely related to the quantum model $`_1`$ in the sense discussed in Section VII. On the other hand, the quantum model $`_3`$ has a well-defined physical interpretation. It belongs to a class of Hamiltonians which finds applications in the field of $`n`$th power squeezed states.

Finally, the operator $`_1`$ , considered in Section VI, has deficiency indices (1,1) and allows a one-parameter family of self-adjoint extensions, each having a purely discrete spectrum on the real line. The spectra of two different extensions have no point in common. Since different self-adjoint extensions correspond to different dynamics, we needed to fix a given dynamics. This has been carried out by choosing the value $`\epsilon =0`$ for the eigenvalue parameter. In this case all the differential equations coming from (69) in all the representations: $`n`$ rep, $`q`$ rep and $`z`$ rep, can be solved exactly. In this case, square-integrable solutions of the eigenvalue equations are explicitly determined.

In order to investigate some properties of the Hamiltonians $`_j`$ ($`j=1,2,3,4),`$ in Section VIII we have introduced the operator $`_5(a^3+a^3),`$ which is connected with $`_3`$, as one can see by using the transformation $`a^{}=ia,`$ $`a^{}=ia^{}.`$ We have built up a set of constants of motion involving $`(_1,_4)`$ and $`(_3,_5).`$ These constants play the role of the arbitrary constants present in the general solution of two equations of the harmonic oscillator-type of frequencies $`1`$ and $`3,`$ respectively, which are satisfied by $`(_1,_4)`$ and $`(_3,_5)`$ (see (93) and (94)).

To conclude, we point out that quadratic (and, more in general, polynomially deformed algebras) take place in quantum optics in relation to the construction of coherent states and in the description of multiphoton processes (see, for examples,  and references therein). Keeping in mind these problems, it should be of interest to deal with the quantization of fouled Hamiltonians (of the generalized oscillator) expressed by polynomials in the canonical variables $`(q,p)`$ of degree higher than three.

Acknowledgments

We would like to express our gratitude to the (anonymous) referee for notable suggestions, especially in regard to the improvement of Section VI. We also thank Dr. G. Esposito, of Istituto Nazionale di Fisica Nucleare, Sezione di Napoli, for helpful discussions.

## XI  Appendix: the Casimir invariant for the quadratic algebra (107)

We remind the reader that a standard form of a quadratic algebra, which is important in the treatment of coherent states of trilinear boson Hamiltonians , is

$$[N_0,N_\pm ]=\pm N_\pm ,[N_+,N_{}]=\pm 2N_0+\delta N_0^2,$$
(112)
where the positive (negative) sign of $`2N_0`$ indicates a polynomially deformed $`su(2)`$ $`(su(1,1)),`$ and $`\delta `$ is a parameter.

The Casimir operator is given by

$$C=N_{}N_++N_0(N_0+1)[1+\frac{\delta }{6}(2N_0+1)].$$
(113)
By setting

$$J_0=a_0N_0+b_0,$$
(114)
$$J_{}=k_1N_{},$$
(115)
$$J_+=k_2N_+,$$
(116)
and choosing, for example,

$$a_0=1,b_0=\frac{i}{2\sqrt{3}},k_1k_2=\frac{i\sqrt{3}}{2},$$
(117)
the commutator relations (107) are converted into

$$[N_0,N_\pm ]=\pm N_\pm ,[N_+,N_{}]=2N_0+\delta N_0^2,$$
(118)
with $`\delta =2i\sqrt{3}.`$ Hence, the Casimir operator (113) reads

$$C=N_{}N_++N_0(N_0+1)[1\frac{i}{\sqrt{3}}(2N_0+1)].$$
(119)