1 Introduction

Variable Gamma-ray Emission Induced by Ultra-High Energy Neutral Beams: Application to 4C +21.35

Abstract

The flat spectrum radio quasar (FSRQ) 4C +21.35 (PKS 1222+216) displays prominent nuclear infrared emission from K dust. A 70 – 400 GeV flare with min variations during half an hour of observations was found by the MAGIC telescopes, and GeV variability was observed on sub-day timescales with the Large Area Telescope on Fermi. We examine 4C +21.35, assuming that it is a source of ultra-high energy cosmic rays (UHECRs). UHECR proton acceleration in the inner jet powers a neutral beam of neutrinos, neutrons and rays from photopion production. The radiative efficiency and production spectra of neutrals formed through photohadronic processes with isotropic external target photons of the broad line region and torus are calculated. Secondary radiations made by this process have a beaming factor , where is the Doppler factor. The pair-production optical depth for rays and the photopion efficiency for UHECR neutrons as they pass through external isotropic radiation fields are calculated. If target photons come from the broad line region and dust torus, large Doppler factors, are required to produce rapidly variable secondary radiation with isotropic luminosity erg s at the pc scale. The -ray spectra from leptonic secondaries are calculated from cascades initiated by the UHECR neutron beam at the pc-scale region and fit to the flaring spectrum of 4C +21.35. Detection of TeV neutrinos from 4C +21.35 or other VHE blazars with IceCube or KM3NeT would confirm this scenario.

Subject headings:
black holes: jets—gamma rays: theory—radiation mechanisms: nonthermal

1. Introduction

Three flat spectrum radio quasars, 3C 279 (Albert et al., 2008; Aleksić et al., 2011a), PKS 1510-089 (Wagner et al., 2010; Cortina et al., 2012), and PKS 1222+216 (Aleksić et al., 2011b), have been detected in the very-high energy (VHE; GeV) regime. This last object, a blazar at redshift also known as 4C +21.35, is remarkable for day-scale GeV flares monitored with the Large Area Telescope (LAT) on the Fermi Gamma-ray Space Telescope in 2010 April and June (Tanaka et al., 2011) that coincided with a period of high infrared activity (Carrasco et al., 2010). Moreover, the Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) telescopes observed a VHE flaring episode during a period of observation lasting about one-half hour on 2010 June 17 that showed factor-of-two variations on a 10 minute time scale (Aleksić et al., 2011b). The rapidly variable rays reached large apparent isotropic -ray luminosities erg s in the VHE band and erg s above 100 MeV.

The intense and highly variable -ray fluxes measured from blazars, including 4C +21.35, demonstrate that they are powerful accelerators of energetic particles. Energetic leptons make intense fluxes of synchrotron emission, with associated rays formed by Compton processes from the same electrons. As shown by Fermi LAT observations (Ackermann et al., 2011), the two-year time-averaged values of of dozens of FSRQs exceed erg s. During flaring episodes, values of erg s are not unusual (the record-setting value is erg s from the November 2010 outburst of 3C 454.3; Abdo et al., 2011a). In addition to relativistic leptons, hadrons are also likely to be accelerated in blazar jets and make secondary radiations by interacting with target photons of the internal and surrounding radiation fields. This emission will contribute to the formation of the spectral energy distribution (SED) (see Böttcher, 2010, for a recent review), and the accelerated hadrons or secondary neutrons could escape from the acceleration region to become ultra-high energy cosmic rays (UHECRs) with energies eV.

Typical quasars have broad optical and ultraviolet emission lines that originate from the broad line region (BLR) near the supermassive black hole. The rapid variability, as short as a few hours for 4C +21.35 as found in the Fermi-LAT data (Foschini et al., 2011a, b), indicates that the high-energy rays are formed close to the black hole within the BLR. But if this is the case, such high-energy photons would be strongly attenuated through the process against scattered line and accretion-disk radiation at the sub-pc scale. We argue that the contradictory behaviors of rapid variability, large luminosity, and a production region distant from the central nucleus can be reconciled if the inner jet accelerates UHECRs that escape as neutrons, which then undergo photo-hadronic breakup in the infrared (IR) field of the dusty torus to make rapidly variable synchrotron rays at GeV – TeV energies on the pc scale. Furthermore, UHE rays made in the inner jet can contribute to the cascade at the pc scale.

The recent observations of 4C +21.35 are summarized in Section 2. An overview of the model, with estimates of the radiative efficiency and secondary photohadronic production spectra, is given in Section 3. A detailed treatment of secondary production of neutrals is presented in Section 4, giving beaming factors and efficiencies of secondaries made in relativistic jets. Production spectra of neutrons and neutrinos by accelerated UHECRs are calculated. In Section 5, formulas for and calculations of pair-production opacity and photohadronic efficiency of rays and neutrons as they travel through the BLR are given. We show at the end of this Section that the BLR is highly opaque to GeV and TeV rays, while TeV photons are still strongly attenuated at the pc scale. Ultra-high energy protons and neutrons escaping from the inner jet are shown to lose significant energy through photopion production at the pc scale. In Section 6, the parsec-scale magnetic field values are derived that allow for sub-10 minute variability of lepton synchrotron radiation from electrons and positrons (hereafter referred to as “electrons”) formed as secondaries of photohadronic neutron production. Numerical fits to the MAGIC data of 4C +21.35 are presented. Discussion and summary are given in Section 7.

We make the simplifying assumption throughout that the radiation fields are locally isotropic. This assumption becomes increasingly poor if the radiation energy density is dominated by anisotropically distributed photons impinging from behind. Thus, provided the emitting region is within the BLR radius () for processes involving atomic-line or reprocessed accretion-disk radiation, or within the pc scale for processes involving the IR radiation field of the dusty torus, these results give an accurate determination of the opacity and photopion efficiency. Note also, that for 4C +21.35 at , the luminosity distance Mpc cm in a CDM cosmology, where the Hubble constant , in units of 100 km s Mpc, is the ratio of matter density to the critical density, and is the ratio of cosmological constant to the critical density.

2. Observations of 4C +21.35

In this section, the relevant observations and radiation fields in the nuclear region of 4C +21.35 are described.

2.1. Observations

The discovery of VHE emission between and 400 GeV with the MAGIC telescopes (Mariotti et al., 2010; Aleksić et al., 2011b) was anticipated by the report of VHE photons found in Fermi LAT data during a high state of 4C +21.35 in April 2010 (Neronov, Semikoz & Vovk, 2010). The source was targeted with the MAGIC telescopes from 3 May 2010. The spectrum of the major VHE flare observed with MAGIC on 17 June 2010, which observed for only minutes, is well described by a single power law with photon number index that is hardened by about one unit when corrected for attenuation with a low-intensity model (Kneiske & Dole, 2010) of the extragalactic background light (EBL) at infrared and optical frequencies. The remarkably short 10 minute variability (in comparison, 3C 279 displayed weak VHE variability on time scales of days; see Albert et al., 2008) points to a production site near the central supermassive black hole where collimation by extreme processes can occur (Tavecchio et al., 2010). Variability associated with the dynamical timescale of a Solar mass black hole is s, where is the Schwarzschild radius of the black hole, whereas variability on time scales shorter than by factors of 10 or more has been observed in the rapidly variable BL Lac objects Mrk 501 (Aharonian et al., 2007), PKS 2155-304 (Aharonian et al., 2007), and Mrk 421 (Fossati et al., 2008). For 4C +21.35, with a black hole mass estimated at (Wang et al., 2004), the corresponding timescale is s, so the VHE measurements indicate a marginally hyper-variable source.

A higher black-hole mass is obtained from recent optical spectroscopy of J1224+2122 by Shaw et al. (2012), using scaling relations derived from reverberation mapping of radio-quiet AGNs. From analysis of the FWHM of the H line, they deduce a black-hole mass for 4C +21.35 of , which does not include uncertainties in applying scaling relations to blazars. This value is considerably higher than the mass reported by Wang et al. (2004), and would relax the energetics based on the Eddington luminosity, but make the short timescale variability more problematic. Here we use the lower mass, and return to this point in Section 7.

The BLR radius can be written in terms of the accretion-disk luminosity as (Ghisellini & Tavecchio, 2008)

(1)

Based on the H luminosity, erg s, Tanaka et al. (2011) derive from the scaling relations of Wang et al. (2004) and Fan et al. (2006) that the BLR luminosity erg s, and that the accretion disk luminosity is an order of magnitude greater, erg s. In comparison, Tavecchio et al. (2011) estimate that erg s, assuming the Swift UVOT spectrum from 4C +21.35 is the optically thick accretion-disk radiation field, in which case cm pc. This represents super-Eddington luminosities, given that erg s unless the black-hole mass is . We use the smaller disk luminosity here, noting that a larger and BLR radius will improve photohadronic efficiency.

More precise relations depend on specific lines and assumptions about the geometry of the BLR (see, e.g., the recent discussion by Foschini, 2011). In particular, we could use the rest-frame 5100 Åor 1350 Åcontinuum, or the high-ionization Lyman line, as seen strongly in, e.g., 3C 454.3 (Bonnoli et al., 2011). For definiteness, we focus on the Ly line to illustrate BLR line effects, which typically carries % of the BLR luminosity (Francis et al., 1991). Approximating the Ly /BLR line luminosity . and defining

(2)

gives a Ly photon field with energy density erg cm.

The magnetic field in the comoving jet frame can be determined if the GeV rays are Compton-scattered external radiation and the lower energy optical radiation is synchrotron radiation made by the same nonthermal electrons making the rays at GeV energies. The Compton dominance parameter is defined as the ratio of the -ray luminosity (assumed to originate from external Compton scattering) to the radio/X-ray synchrotron luminosity. For 4C +21.35, reaches values as large as (see multiwavelength SED using contemporaneous data in Tavecchio et al., 2011). If the external radiation field scattered by these electrons are the Ly photons, then , with , and (e.g., Sikora et al., 2009). Thus the comoving magnetic field for a blazar jet with bulk Lorentz factor is

(3)

Modeling 5 – 35 Spitzer, SDSS, 2MASS, and Swift UVOT data of 4C +21.35, Malmrose et al. (2011) decompose its spectrum into a nonthermal power-law and two-temperature dust model. Hot dust with K radiates erg s from a pc-scale region, and a second warm dust component radiates erg s at K on the same scale. The IR component, which has been proposed as a target in external Compton scenarios (Błażejowski et al., 2000; Sikora et al., 2009), is especially important for UHECR scenarios, noting that the relative energy densities of the BLR, the dust fields, and the CMBR (at ) are, in units of erg cm, , erg s, and , respectively. We use the same parameters as derived by Malmrose et al. (2011) in our model.

2.2. Target radiation fields

Candidate target radiation fields are the cosmic microwave background radiation (CMBR), with dimensionless temperature , 1200 K and 660 K dust ( and , respectively), quasi-thermal scattered 10 eV disk emission (). For Ly line radiation—neglecting broadening, which is unimportant for these calculations— . Using the GZK energy of eV as a yardstick where isobar excitation by 2.725 K CMBR photons is a source of opacity and energy loss, the IR dust and accretion disk/Ly radiation are effective targets for photopion production from cosmic-ray proton and neutron (nucleon) interactions with eV and eV, respectively. The energy densities of the target radiation fields are defined by the usual relation where is the distance from the black hole.4

The luminosity of photopion secondaries is proportional to the effective photon number density for photopion production. For monochromatic line sources, , and for thermal blackbody or graybody radiation fields, , as shown below. Here, is the mean energy of photons in the unit of . The effective photon number density for the CMBR is

(4)

assuming the temperature of the CMBR at the present epoch is 2.725 K. The effective photon density for scattering of the Ly field is

(5)

The BLR radiation field due to the accretion-disk radiation scattered by BLR gas is, when approximated as a graybody, given by

(6)

assuming an effective Thomson scattering depth of . For the two dust radiation fields, the effective photon densities for photopion p interactions are

(7)

and

(8)

Here the luminosity of the IR radiation, assumed to be radiated on a size scale of pc, is denoted erg s. To order of magnitude, cm, and the CMBR is negligible.

For sufficiently energetic cosmic-ray nucleons, the warm dust radiation field provides the densest and most important target photon field for photopion production. For dissipation of energy on the pc scale, only the highest energy nucleons with eV lose energy effectively by photohadronic processes with IR dust photons, because the BLR radiation does not extend to the pc scale, and there is no other sufficiently strong radiation field at the pc scale to extract energy efficiently from these lower energy protons.

3. Model and Estimates

The model is outlined, followed by a simple estimate of photohadronic efficiency and a simple derivation of secondary neutral production spectra made by photohadronic interactions of cosmic-ray jet protons in the inner jet.

3.1. Model description

This model is motivated by the hypothesis that that blazars and radio galaxies accelerate UHECRs (e.g., Mannheim & Biermann, 1992; Berezinsky et al., 2006; Dermer et al., 2009; Murase & Takami, 2009); confirmation of this hypothesis is the detection of UHE neutrinos during blazar flares. Electromagnetic signatures can also provide crucial evidence in support of UHECR production in blazars, and the presence of distinct hadronic emission signatures in the variable -ray SEDs is consistent with the SED shape but not unambiguous (see Böttcher (2010) or, e.g., Abdo et al., 2011b, for Mrk 421). Detection of a multi-TeV radiation signature in extreme high-synchrotron-peaked blazars like 1ES 0229+200 (Murase et al., 2012) at energies where the EBL should suppress this emission would support emission generated by UHECR protons accelerated by blazars (Essey & Kusenko, 2010; Essey et al., 2010, 2011).

Figure 1.— Cartoon of the system. Ultra-high energy neutrals, made through photohadronic production of cosmic-ray protons accelerated in the inner jet of 4C +21.35, escape to form a beam of neutrons, neutrinos, and rays. Subsequent photohadronic and interaction of the neutrons and UHE rays with IR photons from the dust torus induce a high-energy lepton beam that quickly radiates secondary VHE photons to make the variable -ray emission in 4C +21.35. UHE neutrons and ions escaping from the system become UHECRs, and the jet radiation is accompanied by an associated flux of UHE neutrinos.

Fig. 1 shows a cartoon illustrating the underlying basis of this model. UHECR protons and ions are accelerated in the outflowing relativistic plasma formed in the inner jet. Interactions with internal synchrotron and external quasi-isotropic radiation field cause ultra-relativistic protons accelerated in the outflowing jet to undergo photopion losses, with the production of escaping neutrons, rays, and neutrinos (Atoyan & Dermer, 2003). Subsequent interactions of the escaping neutrals with IR torus photons on the pc scale make highly beamed leptons produced as secondaries of interactions of UHE photons, and as secondaries of photopion interactions of UHE neutrons. The relativistic leptons make beamed synchrotron and Compton scattered radiation. These emissions, being highly beamed, will vary on the production timescale determined by the behavior in the inner jet. Off-axis emission from misdirected leptons will not smear and enhance the timescale because it is produced by ultra-relativistic electrons that are not significantly deflected, so the contribution to the observer originates only from the particles traveling directly toward the observer.

3.2. Radiation efficiency estimates

We make simple estimates of the radiative efficiency for the production of photohadronic secondaries when UHE protons interact with photons from isotropic radiation fields external to the jet. In our estimates, the energy density of the external radiation field is denoted and the characteristic spatial size scale of the external radiation field is . For thermal radiation fields, .

Suppose that a relativistic jet moving with bulk Lorentz factor is filled with cosmic-ray protons and ions. The comoving dynamical time for traveling through an external radiation field of is . In the jet frame, the energy density of external radiation field is , and the mean photon energy is . The threshold Lorentz factor for photomeson production by protons with comoving Lorentz factor interacting with photons with energy is given by , where and MeV is the pion mass.

The energy-loss rate of ultra-relativistic protons through photopion processes is therefore

(9)

where b is the product of the inelasticity and the cross section for photomeson losses (Begelman et al., 1990; Atoyan & Dermer, 2003). Note that the energy-loss timescale is proportional to the effective photon density . Employing the approximation of Atoyan & Dermer (2003) here and in the following, we take . The photopion production efficiency is therefore simply

(10)

provided . The Lorentz factor of protons measured in the stationary frame of the black-hole jet system (as if they had escaped the jet) is , so the threshold Lorentz factor is .

Following Equations (5) – (8), we write . For interactions in the BLR with Ly photons or scattered accretion disk radiation, for escaping protons with energy eV. Here, cm. For interactions with IR photons from the dust torus, for protons with eV. For larger dissipation radii, one has , which means that UHECR protons in the jet plasma can dissipate essentially all of their energy into secondaries while traveling through the radiation field of the dust torus. But note that in this case, the UHE proton beam has to be maintained throughout the entire pc-scale length of the external radiation field, in which case the emission would vary on timescales no shorter than .

Thus we see that the efficiency to produce photomeson secondaries through interactions of UHECR jet protons with photons of an external radiation field is several % in the inner jet, counting both line and scattered photons, and noting that erg s (Tanaka et al., 2011), so that cm, from Equation (1). This holds for eV protons. Under optimistic conditions, the photopion efficiency of cosmic-ray protons in a blazar jet is of order unity in interactions with the IR dust photons, but this applies only to protons with escaping energies eV. Note furthermore that this is a conservative estimate of photohadronic efficiency, because it does not take into account internal synchrotron photons, or photons from a sheath-spine structure or from different portions of the jet (compare Ghisellini et al., 2005; Georganopoulos & Kazanas, 2003, for TeV blazars, though similarly structured jets could be found in FSRQs) that could significantly enhance the target radiation field and the photohadronic efficiency.

It is interesting that all the factor dependencies in Equation (10) drop out. This shows that we can equivalently calculate the photohadronic efficiency for a proton traversing the volume of the external radiation field without transforming to and from the comoving frame. We can furthermore make a simple derivation of the secondary neutron production by only considering interactions in the frame of the black hole and external radiation fields.

With a photopion energy-loss cross section of b, the column density of photons above threshold required for unity radiation efficiency is cm. The column density of Ly photons in the BLR is cm. The column density of IR photons from the pc-scale warm dust torus is cm. The IR photon field of the warm dust clearly presents the largest column density for photohadronic production, for the most energetic nucleons above the photopion reaction threshold.

3.3. Simple Derivation of Secondary Neutron Production Spectra

The estimate above implicitly assumes that when protons undergo photopion interactions, they remain protons. In about 50% of the time, however, neutrons are formed which escape from the jetted plasma in a collimated outflow. The escaping neutrons deposit energy throughout the length of the jet by subsequent photohadronic interactions. We now estimate the energy loss rate of protons into neutrons via photopion production, and use that result to make a simple derivation of the secondary neutron production spectrum, anticipating the more detailed derivation in the next section.

The inclusive energy-loss rate of protons into neutrons only is given by

(11)

where is the cross-section ratio, is the fractional energy of the incident proton that is deposited into secondary neutrons, and the neutron multiplicity. For neutron production near threshold, , , , and the maximum photopion cross section is b. Here the proton energy , and .

The rate for a primary proton to lose energy into secondary neutrons is therefore simply

(12)

where is the invariant interaction energy. Solving for a monochromatic radiation source gives

(13)

where .5

The inclusive energy loss by primary protons with Lorentz factor that is transformed into secondary neutrons with Lorentz factor is conserved; therefore

(14)

and . The neutron production luminosity is . For a stationary-frame distribution in the effective volume (where the retarded time is considered), of protons differential in and direction , the secondary production spectrum of neutrons, assuming that the neutrons have the same direction as the original protons, is

(15)

Using the method of Georganopoulos et al. (2001, 2004), we transform the comoving jet frame proton spectrum to the stationary frame jet proton spectrum through the relation

(16)

where the last relation is a consequence of the assumption of proton isotropy in the proper fluid frame, noting that . Thus

(17)

Note the beaming factor arising from the transformation of the proton spectrum, equation (16), and one power each from energy and time.

From Equation (17), we can determine the absolute amount of energy in protons, that, while traveling through a radiation field characterized by , will produce an apparent isotropic luminosity in neutrons erg s. If the spectrum of protons is , that is having equal energy per decade in protons, and extends well above the threshold for neutron production, then the fractional energy is proportional to the ratio of two logarithmic factors, which has a value of order , implying

(18)

There is no fundamental reason why the engine should be Eddington-limited during outbursts, but this provides a fiducial amount of energy that could be generated by the central engine. Therefore we assume that the absolute nonthermal proton energy erg, for the 4C +21.35 flare ( s is the variability timescale of the flare and ), which implies that

(19)

so in order to make an apparent luminosity erg s in neutrons.

The factor limit can be relaxed to if we happen to be viewing exactly down the jet () or if the energy content in protons is not limited by the Eddington luminosity during this outburst, but is instead one or two orders of magnitude larger. Alternately, small is possible if the photomeson production occurs mainly by target photons produced in inner jets (see below). In principle, there is no real reason that such large values of the Doppler factor or factor are not allowed. Arguments based on opacity applied to the July/August 2006 flares from PKS 2155-304 (Aharonian et al., 2007) require (Begelman et al., 2008), while a full synchrotron/SSC model fit, including and EBL effects, implies for PKS 2155-304 and for Mrk 501 with a variability timescale of s (Finke et al., 2008). One-zone leptonic models for the VHE emission from 3C 279 (Albert et al., 2008; Aleksić et al., 2011a) furthermore require large, , bulk Lorentz factors (Böttcher et al., 2009). Standard arguments relating variability time and location of emission site through the expression imply for emission produced at the pc scale, which already undermines naive expectations about and the size of the emitting region. Measurements of superluminal motion never reach such values, but this is based on radio observations, which may be tracking the outflow speed of a radio-emitting sheath rather than a more rapidly moving spine. The necessity of a radical departure from standard models to explain the 4C +21.35 result seems unavoidable in all models (Tavecchio et al., 2011; Nalewajko et al., 2012; Tavecchio et al., 2012); here we construct a model within a framework that is testable by neutrino telescopes.

4. Spectral Powers and Beaming Factors for Photohadronic Secondaries

The system is now treated more carefully. A standard blazar jet scenario is considered (Atoyan & Dermer, 2001). In the comoving frame defined by the existence of a quasi-isotropic particle and randomized magnetic-field distribution moving with Lorentz factor at some angle with respect to the observer, secondaries are formed as a result of photohadronic processes, e.g., neutrons (n), neutrinos (), electrons and positrons (e), photons (), protons (p), and -decay electrons (e) and neutrinos (). The charged secondaries may be trapped in the jet by its magnetic field, but the neutrals (n, , and rays) escape, forming a neutral beam with a beaming factor reflecting the Doppler factor of the jet and the production kinematics (Atoyan & Dermer, 2003).

We start with Equation (2.44) of Dermer & Menon (2009), but instead of treating electron-photon scattering, here we treat secondaries of protons with Lorentz factor made in photopion interactions with photons with dimensionless energy . The technique is, however, general and can be adapted to other secondaries, noting that the dimensionless secondary energy in units will be written as a fraction of of the primary proton energy . Thus for secondary photons, leptons or neutrinos, . For secondary neutrons, we can directly write, because , the neutron production spectrum in terms of the secondary neutron Lorentz factor . The direction-dependent emissivity for the production of neutrons through the process becomes

(20)

specialized to ultra-relativistic protons, . Here is the angle between the directions of the interacting proton and photon. The approximation we use for the differential cross section for photo-pion processes is a sum of step functions, given by

(21)

Here the relativistic invariant is , , and and are the cosine angle and azimuthal angle, respectively, of the photon or (with subscript “p”) proton. The term is the fractional cross section written as a ratio of the value of the step-function cross section to the maximum cross section b. The effective multiplicity for segment between invariant energies is denoted , and is the mean energy fraction of a secondary neutron compared to the original proton energy for interactions associated with segment .

Substituting Equation (21) into (20) and solving gives

(22)

where

and , taking without loss of generality.

Using Equation (16), we have

(23)

where .

We now specialize to the case of a surrounding external radiation field that is isotropic in the stationary frame of the black-hole/accretion-disk system. Therefore , and we can let without loss of generality. After some straightforward manipulations, we obtain

(24)

where , and

Or, we have

(25)

4.1. External isotropic monochromatic radiation field

For a monochromatic external radiation field, we write . In this case,

(26)

This can also be written as

(27)

The single step-function approximation with takes the form

(28)

with , , and . This confirms the approximate derivation, Equation (17), and also shows that secondary production of photohadronic secondaries has a beaming factor .

4.2. External graybody radiation field

Figure 2.— Functions (red; Equation (33)) and (blue; Equation (34)) with asymptotes (dashed and dotted lines). Inset shows functions on a log scale.

The graybody radiation field, defined as a blackbody radiation field times the graybody factor giving the ratio of energy densities of the radiation field under consideration to that of a blackbody radiation field with the same effective temperature is expressed as a spectral energy density in the form

(29)

The electron Compton wavelength cm. The energy density of blackbody radiation is

(30)

and , by definition.

Substituting Equation (29) into (24) and solving gives

(31)

where

(32)
(33)

and