text:
# A spin-cell for spin current

## Abstract

We propose and investigate a spin-cell device which provides the necessary spin-motive force to drive a spin current for future spintronic circuits. Our spin-cell have four basic characteristics: (i) it has two poles so that a spin current flows in from one pole and out from the other pole, this way a complete spin-circuit can be established; (ii) it has a source of energy to drive the spin current; (iii) it maintains spin coherence so that a sizable spin current can be delivered; (iv) it drives a spin current without a charge current. The proposed spin-cell for spin current should be realizable using technologies presently available.

Traditional electronics is based on the flow of charge: the spin of electron is ignored. The emerging technology of spintronics will make the leap such that the flow of spin—in addition to charge, will be used for electronic applications. A spin current is produced by the motion of spin-polarized electrons, therefore spin current is typically associated with spin-polarized charge current. Nevertheless, if one can generate an ideal situation as shown in Fig.(1a) where spin-up electrons move to the right while an equal number of spin-down electrons move to the left, then there will be no net charge current because $`I_ee(I_{}+I_{})=0`$ where $`eI_{},eI_{}`$ are charge currents due to spin-up and spin-down electrons, respectively. There will be, however, a finite spin current: $`I_s\frac{\mathrm{}}{2}(I_{}I_{})`$ where $`\mathrm{}`$ is the reduced Planck constant. Considering the interesting and important future perspective of spin-current circuit, it is crucial to have a spin-cell that satisfies the four characteristics discussed in the Abstract and it produces the flow pattern of Fig.(1a). In this paper we theoretically propose and analyze such a spin-cell.

Our spin-cell is schematically shown in Fig.(1b). It consists of a double quantum-dot (QD) fabricated in two-dimensional electron gas (2DEG) with split gate technology, and each QD is contacted by an electrode. Note that no magnetic material is involved. The two QDs and their associated contacts to the electrodes serve as the “positive/negative” poles of the spin-cell. The two electrodes maintain the same electrochemical potential $`\mu _L=\mu _R`$ (i.e. no bias voltage is applied on them). The size of the spin-cell structure is assumed to be within the spin coherence length which can be as long as many microns for 2DEG. We control the QD energy levels by gate voltages $`V_{g\alpha }`$ where $`\alpha =L,R`$ indicates the left/right QD. Both QD levels are controlled by an overall gate voltage $`V_g`$, see gate arrangements in Fig.(1b). In order to distinguish spin-up electrons from spin-down electrons, a spatially non-uniform external magnetic field $`B_\alpha `$ is applied to the two QDs–perpendicular to the QD plane. An extreme case of non-uniformity is $`B_R=B_L`$, i.e. equal in value but opposite in direction. This particular magnetic field distribution is not necessary at all for the operation of our spin-cell, but it helps us to discuss its physics. Finally, the energy source of our spin-cell is provided by shining a microwave radiation with strengths $`\mathrm{\Delta }_L/\mathrm{\Delta }_R`$ for the left/right QDs. Because, typically, the microwave frequency is far less than the plasma frequency of the material covering the QDs, the effect of the microwave field is to induce a high frequency potential variation $`\mathrm{\Delta }_{L/R}\mathrm{cos}\omega t`$ in the left/right QD and their leads. When $`\mathrm{\Delta }_L\mathrm{\Delta }_R`$, a time-dependent potential difference, $`\mathrm{\Delta }\mathrm{cos}\omega t=(\mathrm{\Delta }_L\mathrm{\Delta }_R)\mathrm{cos}\omega t`$, exists between the two QDs. An a.c. electric field $`𝐄(t)`$ in the middle barrier is therefore established due to the microwave radiation (see Fig.(1c)). Then, electrons can absorb photons when they pass the middle barrier of the device. The establishment of $`𝐄(t)`$ across the two QDs is necessary for our spin-cell to work, here we use a non-uniform microwave radiation to achieve this effect as has already been carried out experimentally, but other possibilities also exist.

Before we present theoretical and numerical results of the device in Fig.(1b), we first discuss why it works as a spin-cell. The physics is summarized in Fig.(1c). To be specific, let $`B_R`$ point to $`z`$ direction and $`B_L`$ to $`+z`$ direction. Due to Zeeman effect, a spin-degenerate level $`ϵ_R`$ on the right QD is now split into spin-down/up levels $`ϵ_R<ϵ_R`$. On the left QD, it is $`ϵ_L<ϵ_L`$. Electrons in the electrodes can now tunnel into the QD: on the right a spin-down electron is easier to tunnel because level $`ϵ_R`$ is lower, while a spin-up electron is easier to tunnel into the left QD. Once levels $`ϵ_R,ϵ_L`$ are occupied, the charging energies $`U_R,U_L`$ of the two QDs push the other two levels $`ϵ_R,ϵ_L`$ to higher energies $`ϵ_R+U_R,ϵ_L+U_L`$, and the energy level positions indicated by the solid horizontal lines of Fig.(1c) are established. Next, the spin-down electron on the right QD can absorb a photon and make a transition to the level at $`ϵ_L+U_L`$ on the left QD: afterwards it easily flows out to the left electrode because $`ϵ_L+U_L>\mu _L`$. This process is indicated as $`A`$. Similarly the spin-up electron on the left QD flows out to the right electrode after absorption of a photon, indicated by $`A+`$. This way, driven by the potential variations of the QD induced by the microwave field, a spin-down electron flows to the left while a spin-up electron flows to the right of the spin-cell, and the continuation of the $`A,A+`$ processes generates a DC spin current that flows from the left electrode, through the spin-cell, and out to the right electrode. Clearly, if the two processes are absolutely equivalent, there will be no charge current and only a spin current. Finally, since the spin-motive force is provided by a time-dependent change of the electronic potential landscape of the QD, there is no spin flip mechanism and the spin current flowing through the spin-cell is conserved, i.e. $`I_{s,L}=I_{s,R}=I_s`$. Our device then satisfies the four characteristics of a spin-cell discussed in the Abstract.

The last paragraph discusses the operation principle of the spin-cell for spin current, but there are other interesting device details which can only be obtained by detailed theoretical and numerical analysis for which we now turn. The spin-cell of Fig.(1b) is described by the following Hamiltonian:

$`H`$ $`={\displaystyle \underset{\alpha \sigma }{}}\left[ϵ_\alpha +W_\alpha (t)(1/2)\sigma g\mu B_\alpha \right]d_{\alpha \sigma }^{}d_{\alpha \sigma }`$ (3)
$`+{\displaystyle \underset{\alpha }{}}U_\alpha d_\alpha ^{}d_\alpha d_\alpha ^{}d_\alpha +{\displaystyle \underset{k\sigma \alpha }{}}\left[ϵ_{\alpha k}+W_\alpha (t)\right]a_{\alpha k\sigma }^{}a_{\alpha k\sigma }`$
$`+{\displaystyle \underset{k\sigma \alpha }{}}[t_{\alpha k}a_{\alpha k\sigma }^{}d_{\alpha \sigma }+H.c.]+{\displaystyle \underset{\sigma }{}}[t_Cd_{L\sigma }^{}d_{R\sigma }+H.c.]`$
where $`a_{\alpha k\sigma }^{}`$($`a_{\alpha k\sigma }`$) and $`d_{\alpha \sigma }^{}`$($`d_{\alpha \sigma }`$) are creation (annihilation) operators in the electrode $`\alpha `$ and the dot $`\alpha `$, respectively. The left and right QD includes a single energy level $`ϵ_\alpha `$, but has spin index $`\sigma `$ and intradot Coulomb interaction $`U_\alpha `$. To account for the magnetic field $`B`$, the left/right QD’s single particle energy has a term $`(1/2)\sigma g\mu B_\alpha `$, in which we have required a different magnetic field strength for the two QDs, i.e. $`B_LB_R`$. $`t_C`$ and $`\mathrm{\Gamma }_\alpha 2\pi _k|t_{\alpha k}|^2\delta (ϵϵ_{\alpha k})`$ describe the coupling strength between the two QDs, and between electrode $`\alpha `$ and its corresponding QD, respectively. The microwave irradiation is given by $`W_\alpha (t)=\mathrm{\Delta }_\alpha \mathrm{cos}\omega t`$ and it produces an adiabatic change for the single particle energy. Here we permit the microwave field to irradiate the entire device including the electrodes, and we require a difference in the radiation strength $`\mathrm{\Delta }_L\mathrm{\Delta }_R`$.

Our theoretical analysis of the spin-cell is based on standard Keldysh nonequilibrium Green’s function (NEGF) theory which we briefly outline here. First, we perform an unitary transformation of the Hamiltonian with an unitary operator $`U(t)=exp\left\{i_0^t𝑑t^{}\underset{\alpha }{}W_\alpha (t^{})\widehat{D}\right\}`$, where $`\widehat{D}\underset{k\sigma }{}a_{\alpha k\sigma }^{}a_{\alpha k\sigma }+\underset{\sigma }{}d_{\alpha \sigma }^{}d_{\alpha \sigma }`$. The Hamiltonian $`H`$ is transformed to the following form,

$`H`$ $`=`$ $`{\displaystyle \underset{\alpha \sigma }{}}\left[ϵ_\alpha \sigma g\mu B_\alpha /2\right]d_{\alpha \sigma }^{}d_{\alpha \sigma }+{\displaystyle \underset{\alpha }{}}U_\alpha d_\alpha ^{}d_\alpha d_\alpha ^{}d_\alpha `$ (6)
$`+{\displaystyle \underset{k\sigma \alpha }{}}ϵ_{\alpha k}a_{\alpha k\sigma }^{}a_{\alpha k\sigma }+{\displaystyle \underset{k\sigma \alpha }{}}[t_{\alpha k}a_{\alpha k\sigma }^{}d_{\alpha \sigma }+H.c.]`$
$`+{\displaystyle \underset{\sigma }{}}[t_Ce^{i_0^t𝑑t^{}\mathrm{\Delta }cos\omega t^{}}d_{L\sigma }^{}d_{R\sigma }+H.c.],`$
where $`\mathrm{\Delta }\mathrm{\Delta }_L\mathrm{\Delta }_R`$. In (6), we take the last term which explicitly depends on time $`t`$ as the interacting part $`H_I`$, and the remaining part as $`H_oHH_I`$. The Green’s function of $`H_o`$, $`g^r(ϵ)`$, can be easily obtained with a decoupling approximation at the Hartree level,

$$g_{\alpha \alpha \sigma }^r(ϵ)=\frac{ϵ_{\alpha \sigma }^{}+U_\alpha n_{\alpha \overline{\sigma }}}{(ϵϵ_{\alpha \sigma })ϵ_{\alpha \sigma }^{}+\frac{i}{2}\mathrm{\Gamma }_\alpha (ϵ_{\alpha \sigma }^{}+U_\alpha n_{\alpha \overline{\sigma }})},$$
(7)
where $`ϵ_{\alpha \sigma }^{}ϵϵ_{\alpha \sigma }U_\alpha `$, $`ϵ_{\alpha \sigma }ϵ_\alpha \sigma g\mu B_\alpha /2`$, and $`n_{\alpha \overline{\sigma }}`$ is the time-averaged intradot electron occupation number at the state $`\overline{\sigma }`$ in the $`\alpha `$-QD which we solve self-consistently. It is worth mentioning that $`g_{\alpha \alpha \sigma }^r(ϵ)`$ in Eq.(7) has two resonances: one is at $`ϵ_{\alpha \sigma }`$ while its associated state at $`ϵ_{\alpha \overline{\sigma }}`$ is empty; the other resonance is at $`ϵ_{\alpha \sigma }+U_\alpha `$ while its associated state $`ϵ_{\alpha \overline{\sigma }}`$ is occupied. Notice, in $`H_o`$ the left part of the spin-cell (i.e. the left-lead and the left-QD) is not coupled with the right part of the spin-cell, therefore they are in equilibrium respectively. Hence the Keldysh Green’s function $`g_{\alpha \alpha \sigma }^<(ϵ)`$ for $`H_o`$ can be solved from the fluctuation-dissipation theorem: $`g_{\alpha \alpha \sigma }^<(ϵ)=f_\alpha \left[g_{\alpha \alpha \sigma }^r(ϵ)g_{\alpha \alpha \sigma }^a(ϵ)\right].`$ With these preparations, the Green’s function $`G^r`$ and $`G^<`$ of the total Hamiltonian $`H`$ can be solved. In particular, we calculate $`G_{\alpha \beta \sigma }^r(t,t^{})i\theta (tt^{})<\{d_{\alpha \sigma }(t),d_{\beta \sigma }^{}(t^{})\}>`$ by iterating the Dyson equation. In Fourier space, the Dyson equation can be reduced to

$$𝐆_{\sigma ;mn}^r(ϵ)=𝐠_{\sigma ;mn}^r(ϵ)+\underset{k}{}𝐆_{\sigma ;mk}^r(ϵ)𝚺_{\sigma ;kn}^r(ϵ)𝐠_{\sigma ;nn}^r(ϵ),$$

where $`𝐆_{\sigma ,mn}^r(ϵ)𝐆_{\sigma ,nm}^r(ϵ+m\omega )`$, and the quantity $`G_n(ϵ)`$ is the Fourier expansion of $`G(t,t^{})`$. The retarded self-energy $`𝚺_{\sigma ;kn}^r(ϵ)`$ is the Fourier transform of $`𝚺_\sigma ^r(t_1,t_2)`$ where $`\mathrm{\Sigma }_{LR\sigma }^r(t_1,t_2)=\mathrm{\Sigma }_{RL\sigma }^r(t_1,t_2)=\delta (t_1t_2)t_Cexp[i_0^{t_1}𝑑t^{}\mathrm{\Delta }cos\omega t^{}]`$, and $`\mathrm{\Sigma }_{LL\sigma }^r=\mathrm{\Sigma }_{RR\sigma }^r=0`$. We obtain $`\mathrm{\Sigma }_{LR\sigma ,mn}^r(ϵ)=\mathrm{\Sigma }_{RL\sigma ,nm}^r(ϵ)=t_CJ_{nm}(\frac{\mathrm{\Delta }}{\omega })`$. The Green’s function $`𝐠_{\sigma ;mn}^r(ϵ)`$ is: $`g_{\alpha \beta \sigma ;mn}^r(ϵ)=\delta _{\alpha \beta }\delta _{mn}g_{\alpha \alpha \sigma }^r(ϵ+m\omega )`$. Then $`𝐆_{\sigma ;mn}^r(ϵ)`$ can be solved from the above Dyson equation:

$$G_{\alpha \alpha \sigma ,mn}^r(ϵ)=\delta _{mn}/\left[\left(g_{\alpha \alpha \sigma ,mm}^r\right)^1A_{\overline{\alpha }\sigma ,mm}\right],$$

$$G_{\alpha \overline{\alpha }\sigma ,mn}^r(ϵ)=G_{\alpha \alpha \sigma ,mm}^r(ϵ)\mathrm{\Sigma }_{\alpha \overline{\alpha }\sigma ,mn}^rg_{\overline{\alpha }\overline{\alpha }\sigma ,nn}^r(ϵ),$$

where $`A_{\alpha \sigma ,mm}(ϵ)\underset{k}{}|t_C|^2J_{km}^2\left(\frac{\mathrm{\Delta }}{\omega }\right)g_{\alpha \alpha \sigma ;kk}^r(ϵ).`$ Afterwards, the total Keldysh Green’s function $`G_{\alpha \beta \sigma }^<(t,t^{})i<d_{\beta \sigma }^{}(t^{})d_{\alpha \sigma }(t)>`$ is easily obtained from the Keldysh equation. Finally, we obtain the time-averaged current in lead $`\alpha `$ from

$`I_{\alpha \sigma }`$ $``$ $`<I_{\alpha \sigma }(t)>=Im{\displaystyle }(dϵ/2\pi )\mathrm{\Gamma }_\alpha (ϵ)[G_{\alpha \alpha \sigma ,00}^<(ϵ)`$ (8)
$`+`$ $`2f_\alpha (ϵ)G_{\alpha \alpha \sigma ,00}^r(ϵ)],`$ (9)
and the self-consistent equation for the intradot occupation number $`n_{\alpha \sigma }`$: $`n_{\alpha \sigma }=i(dϵ/2\pi )G_{\alpha \alpha \sigma ,00}^<(ϵ)`$.

Fig.(2) shows the calculated charge current $`I_e`$ (in units of $`e`$) and the spin current $`I_s`$ (units of $`\mathrm{}/2`$) versus the gate voltage $`V_{gR}`$ at different microwave frequency $`\omega `$. $`I_e`$ shows a positive peak due to the $`A+`$ process and a negative peak by the $`A`$ process (see Fig.(1c)), but $`I_s`$ has two positive peaks. As we tune the gate voltage $`V_{gR}`$, the right QD level is shifted so that when $`\mathrm{}\omega =ϵ_L+U_Lϵ_R`$, the $`A`$ process occurs with high probability leading to a positive peak in $`I_s`$ and a negative peak in $`I_e`$. On the other hand we get positive peaks in both $`I_e`$ and $`I_s`$ when $`\mathrm{}\omega =ϵ_R+U_Rϵ_L`$, for the the $`A+`$ process.

The peak positions in $`I_e,I_s`$ due to the $`A\pm `$ processes shift linearly with the microwave frequency $`\omega `$, as shown by the dotted lines in Fig.(2). Eventually, at a special frequency indicated by $`A`$, i.e. when $`\mathrm{}\omega ^{}=ϵ_R+U_Rϵ_L=ϵ_L+U_Lϵ_R`$, the two peaks overlap so that the net charge current $`I_e`$ cancels exactly due to the cancellation of the $`A\pm `$ processes, at the same time the spin current $`I_s`$ doubles its value. At this special frequency, the full operation of the spin-cell occurs so that a spin current is driven across the spin-cell, from left electrode to the right electrode, without a charge current. If we connect the spin-cell to complete an external circuit, a spin current will be driven and continue to flow across the spin-cell into the circuit. On the other hand, if we let the two poles of the spin-cell open, although $`I_s`$ must be zero, a spin-motive force in the two poles of the spin cell will still be induced so that chemical potential $`\mu _\alpha \mu _\alpha `$. For example, in the case of Fig.(1c), an open circuit will lead to $`\mu _L<\mu _L`$ and $`\mu _R>\mu _R`$.

In the following we focus on the spin-cell operation by fixing gate voltage $`V_{gR}=0.45`$ which is its value at point A of Fig.(2). We investigate $`I_e,I_s`$ as functions of the overall gate potential $`V_g`$ (Fig.(3a)), magnetic field $`g\mu B_L`$ (Fig.(3b)), and frequency $`\omega `$ (Fig.(3c)). Different curves in Fig.(3) correspond to different microwave strength $`\mathrm{\Delta }\mathrm{\Delta }_L\mathrm{\Delta }_R`$. In all situations $`I_e0`$ and we do not discuss it anymore. Fig.(3c) shows that $`I_s`$ has several peaks and dips when we vary $`\omega `$: the large peak indicated by $`A`$ is the spin-cell operation discussed above, but peaks at $`C`$ and $`D`$ correspond to double- and triple-photon processes which connect the $`A\pm `$ transitions of Fig.(1c). The dip at $`B`$ originates from less probable transitions connecting levels indicated by the dashed lines of Fig.(1c), while the dip at $`E`$ is its two-photon process. Now, fixing $`\omega `$ at $`\omega ^{}`$, i.e. at the spin-cell operation point A, the value of $`I_s`$ can be tuned by the overall gate voltage $`V_g`$ as shown in Fig.(3a). However, $`I_s`$ keeps large values for a wide range of $`V_g`$: this range is in the Coulomb interaction scale $`U/e`$! This is important, because in an experimental situation any background charge or environmental effect near the spin-cell may alter the overall potential, and Fig.(3a) shows that the spin-cell operation is not critically altered by this effect. When $`V_g`$ becomes very large so that $`ϵ_L+U_L`$ and $`ϵ_R+U_R`$ is below the chemical potential $`\mu `$, or $`ϵ_L`$ and $`ϵ_R`$ is above $`\mu `$, $`I_s`$ diminishes because the $`A\pm `$ processes can no longer occur (see Fig.(1c)). Finally, a very important result is shown in Fig.(3b), where we fixed $`g\mu B_R=0.4`$ while varying $`g\mu B_L`$ at the spin-cell operation point $`A`$. Fig.(3b) shows clearly that $`I_s`$ increases with an increasing difference of ($`B_LB_R`$): $`I_s=0`$ identically when $`B_L=B_R`$ if $`U_L=U_R`$, or $`I_s0`$ if $`U_LU_R`$. However, Fig.(3b) demonstrates that we only need a slight difference in $`B_L`$ and $`B_R`$, at a scale of the coupling constant $`\mathrm{\Gamma }_\alpha `$, to generate a substantial $`I_s`$. The most important fact is that $`B_L`$ and $`B_R`$ do not have to point to opposite directions which is experimentally difficult to do. In fact, if the two QD’s are fabricated with different materials so that the $`g`$factors are different, one can actually use an uniform magnetic field throughout.

The proposed spin-cell for spin current should be experimentally feasible using present technologies. First, the double-QD structures can and have been fabricated by several laboratories. Second, microwave assisted quantum transport measurements have recently been reported. In particular, the asymmetrical microwave radiation on double-QD device (i.e. $`\mathrm{\Delta }_L\mathrm{\Delta }_R`$) has already been carried out experimentally. Third, the asymmetric magnetic field should be feasible as we have discussed above. If one takes $`f=\omega /2\pi =50`$GHz, arranges the corresponding $`U`$($`\mathrm{}\omega `$)$`0.2meV`$, and fixes the temperature scale $`K_BT`$ and coupling $`\mathrm{\Gamma }_\alpha `$ to be twenty times less than $`U`$ as in typical QD experiments, ı.e. $`k_BT=100mK`$ and $`\mathrm{\Gamma }=10\mu eV`$, the corresponding magnetic field difference is ($`g\mu (B_LB_R)\mathrm{\Gamma }`$) $`|B_LB_R|0.16/g`$ tesla. These QD parameters have already been realized by present technology. Finally, it is not difficult to show that by adjusting the gate voltages one can easily calibrate the spin-cell operating point.

Acknowledgments: We gratefully acknowledge financial support from NSERC of Canada, FCAR of Quebec (Q.S., H.G), and a RGC grant from the SAR Government of Hong Kong under grant number HKU 7091/01P (J.W.). H.G. thanks Dr. Junren Shi for a discussion on photon assisted tunneling. We gratefully acknowledge Dr. Baigeng Wang for many discussions and inputs on the physics of spin current.

Appendix A

In this appendix, we give a detailed derivation process of the unitary transformation in which the Hamiltonian (1) is transformed to Hamiltonian (2) in the manuscript. We take the unitary matrix $`U(t)`$ as:

$$U(t)=exp\left\{i_0^t𝑑t^{}\underset{\alpha }{}\left[W_\alpha (t^{})\left(\underset{k\sigma }{}a_{\alpha k\sigma }^{}a_{\alpha k\sigma }+\underset{\sigma }{}d_{\alpha \sigma }^{}d_{\alpha \sigma }\right)\right]\right\}$$
$`(A1)`$
Under this unitary transformation, the operator $`X_\alpha `$ ($`X_\alpha `$ represents $`a_{\alpha k\sigma }`$ and $`d_{\alpha \sigma }`$ ) and $`X_\alpha ^{}`$ transform into:

$$U(t)X_\alpha U^{}(t)=X_\alpha exp\left[i_0^t𝑑t^{}W_\alpha (t^{})\right]$$
$`(A2)`$
$$U(t)X_\alpha ^{}U^{}(t)=X_\alpha ^{}exp\left[i_0^t𝑑t^{}W_\alpha (t^{})\right]$$
$`(A3)`$
Because unitary matrix $`U(t)`$ is dependent of time $`t`$, the Hamiltonian $`H`$ after the unitary transformation is:

$$H_{new}(t)=U(t)H(t)U^{}(t)U(t)i\frac{}{t}U^{}(t)$$
$`(A4)`$
Using Eqs.(A2) and (A3), we easily obtain the new Hamiltonian (2) as given in the manuscript. We emphasize that this transform is done exactly and no approximation is involved.

Appendix B

In this Appendix, we present a detailed derivation for the retarded Green’s function $`g_{\alpha \alpha \sigma }^r(ϵ)`$ (Eq.(3) in the manuscript) of the Hamiltonian $`H_0`$. Here, we solve $`g_{\alpha \alpha \sigma }^r(ϵ)`$ by using the equation of motion method:

$$ϵ<<\widehat{A}|\widehat{B}>>^r=<\{\widehat{A},\widehat{B}\}>+<<[\widehat{A},H_0]|\widehat{B}>>^r$$
$`(B1)`$
where $`\widehat{A}`$ and $`\widehat{B}`$ represent arbitrary operators, and $`<<\widehat{A}|\widehat{B}>>^r`$ is the Fourier expansion of the retarded Green’s function $`i\theta (t)<\{\widehat{A}(t),\widehat{B}(0)\}>`$. Because $`H_0`$ is independent of time $`t`$, $`<<\widehat{A}(t)|\widehat{B}(t^{})>>^r`$ only depends on the time difference $`tt^{}`$. Appling equation of motion to the Green function $`g_{\alpha \alpha \sigma }^r=<<d_{\alpha \sigma }|d_{\alpha \sigma }^{}>>^r`$, we have:

$$(ϵϵ_{\alpha \sigma })<<d_{\alpha \sigma }|d_{\alpha \sigma }^{}>>^r=1+\underset{k}{}t_{k\alpha }^{}<<a_{\alpha k\sigma }|d_{\alpha \sigma }^{}>>^r+U_\alpha <<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r$$
$`(B2)`$
and

$$<<a_{\alpha k\sigma }|d_{\alpha \sigma }^{}>>^r=\frac{t_{k\alpha }}{ϵϵ_{k\sigma }+i0^+}<<d_{\alpha \sigma }|d_{\alpha \sigma }^{}>>^r$$
$`(B3)`$
Substitute Eq.(B3) into Eq.(B2), we obtain:

$$(ϵϵ_{\alpha \sigma }\mathrm{\Sigma }_\alpha ^r)<<d_{\alpha \sigma }|d_{\alpha \sigma }^{}>>^r=1+U_\alpha <<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r$$
$`(B4)`$
where $`\mathrm{\Sigma }_\alpha ^r\underset{k}{}\frac{|t_{k\alpha }|^2}{ϵϵ_{k\sigma }+i0^+}=\frac{i}{2}\mathrm{\Gamma }_\alpha `$ is the self-energy function.

At this point, if one takes a decoupling approximation for the new Green’s function $`<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r`$ in Eq.(B4), i.e.,

$$<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r=n_{\alpha \overline{\sigma }}<<d_{\alpha \sigma }|d_{\alpha \sigma }^{}>>^r$$
$`(B5)`$
This is the mean field approximation.

Rather, we continue apply the equation of motion to $`<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r`$, and we have:

$$(ϵϵ_{\alpha \sigma }U)<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r=n_{\alpha \overline{\sigma }}$$

$$+\underset{k}{}t_{\alpha k}^{}<<a_{\alpha k\sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r$$

$$\underset{k}{}t_{\alpha k}<<d_{\alpha \sigma }a_{\alpha k\overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r$$

$$+\underset{k}{}t_{\alpha k}^{}<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}a_{\alpha k\overline{\sigma }}|d_{\alpha \sigma }^{}>>^r$$
$`(B6)`$
In this equation, three new Green’s functions $`<<a_{\alpha k\sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r`$, $`<<d_{\alpha \sigma }a_{\alpha k\overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r`$, and $`<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}a_{\alpha k\overline{\sigma }}|d_{\alpha \sigma }^{}>>^r`$, emerge. We now take the decoupling approximation:

$$<<a_{\alpha k\sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r=n_{\alpha \overline{\sigma }}<<a_{\alpha k\sigma }|d_{\alpha \sigma }^{}>>^r$$

$$<<d_{\alpha \sigma }a_{\alpha k\overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r=<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}a_{\alpha k\overline{\sigma }}|d_{\alpha \sigma }^{}>>^r=0$$
$`(B7)`$
Eq.(B6) reduces into:

$$(ϵϵ_{\alpha \sigma }U)<<d_{\alpha \sigma }d_{\alpha \overline{\sigma }}^{}d_{\alpha \overline{\sigma }}|d_{\alpha \sigma }^{}>>^r=n_{\alpha \overline{\sigma }}+\mathrm{\Sigma }_\alpha n_{\alpha \overline{\sigma }}<<d_{\alpha \sigma }|d_{\alpha \sigma }^{}>>^r$$
$`(B8)`$
Then from Eqs.(B4) and (B8), we can solve the Green’s function $`g_{\alpha \alpha \sigma }^r(ϵ)=<<d_{\alpha \sigma }|d_{\alpha \sigma }^{}>>^r`$ as:

$$g_{\alpha \alpha \sigma }^r(ϵ)=\frac{ϵ_{\alpha \sigma }^{}+U_\alpha n_{\alpha \overline{\sigma }}}{(ϵϵ_{\alpha \sigma })ϵ_{\alpha \sigma }^{}\mathrm{\Sigma }_\alpha ^r(ϵ_{\alpha \sigma }^{}+U_\alpha n_{\alpha \overline{\sigma }})}$$
$`(B9)`$
This is the Eq.(3) in the manuscript.

If one applies the equation of motion one more time, on the three new Green’s functions of Eq.(B6), then takes a decoupling approximation, the Green function is obtained to be:

$$g_{\alpha \alpha \sigma }^r(ϵ)=\frac{1+U_\alpha A_{\alpha \sigma }n_{\alpha \overline{\sigma }}}{ϵϵ_{\alpha \sigma }+\frac{i}{2}\mathrm{\Gamma }_\alpha +U_\alpha A_{\alpha \sigma }(\mathrm{\Sigma }_{\alpha \overline{\sigma }}^{(a)}+\mathrm{\Sigma }_{\alpha \overline{\sigma }}^{(b)})}$$
$`(B10)`$
where the definition of $`A_{\alpha \sigma }`$, $`\mathrm{\Sigma }_{\alpha \overline{\sigma }}^{(a)}`$ and $`\mathrm{\Sigma }_{\alpha \overline{\sigma }}^{(b)}`$, are given in Phys. Rev. B 64,153306 (2001). This solution has a Kondo resonance and it has been applied in some previous work for investigating the Kondo effect.

Appendix C

In this Appendix, we give a detailed derivation for the Green’s function $`G_{\alpha \beta \sigma ,mn}^r(ϵ)`$ of the Hamiltonian $`H`$. First, from the Dyson equation, we have:

$$G_{LL\sigma ,mn}^r=g_{LL\sigma ,mn}^r+\underset{k}{}G_{LR\sigma ,mk}^r\mathrm{\Sigma }_{RL\sigma ,kn}^rg_{LL\sigma ,nn}^r$$
$`(C1)`$
$$G_{LR\sigma ,mn}^r=\underset{k}{}G_{LL\sigma ,mk}^r\mathrm{\Sigma }_{LR\sigma ,kn}^rg_{RR\sigma ,nn}^r$$
$`(C2)`$
where we have suppressed the argument $`ϵ`$. Substitute Eq.(C2) into Eq.(C1), we obtain:

$$G_{LL\sigma ,mn}^r=g_{LL\sigma ,mn}^r+\underset{k_1,k_2}{}G_{LL\sigma ,mk_1}^r\mathrm{\Sigma }_{LR\sigma ,k_1k_2}^rg_{RR\sigma ,k_2k_2}^r\mathrm{\Sigma }_{RL\sigma ,k_2n}^rg_{LL\sigma ,nn}^r$$
$`(C3)`$
Introducing notation $`A_{R\sigma ,mn}(ϵ)`$,

$$A_{R\sigma ,mn}(ϵ)\underset{k}{}\mathrm{\Sigma }_{LR\sigma ,mk}^rg_{RR\sigma ,kk}^r\mathrm{\Sigma }_{RL\sigma ,kn}^r$$
$`(C4)`$
then the Eq.(C3) reduces to:

$$G_{LL\sigma ,mn}^r=g_{LL\sigma ,mn}^r+\underset{k_1}{}G_{LL\sigma ,mk_1}^rA_{R\sigma ,k_1n}g_{LL\sigma ,nn}^r$$

$$=g_{LL\sigma ,mn}^r+\underset{k_1}{}g_{LL\sigma ,mk_1}^rA_{R\sigma ,k_1n}g_{LL\sigma ,nn}^r$$

$$+\underset{k_1k_2}{}g_{LL\sigma ,mk_1}^rA_{R\sigma ,k_1k_2}g_{LL\sigma ,k_2k_2}^rA_{R\sigma ,k_2n}g_{LL\sigma ,nn}^r+\mathrm{}$$
$`(C5)`$
In the following we take an approximation as our previous work \[Eq.(29) in Phys. Rev B 59,13126 (1999)\]:

$$A_{R\sigma ,k_1n}g_{LL\sigma ,k_1k_1}^rg_{LL\sigma ,nn}^r=A_{R\sigma ,k_1n}\delta _{k_1n}\left[g_{LL\sigma ,nn}^r\right]^2$$
$`(C6)`$
This approximation is reasonable when $`\mathrm{}\omega >>max(\mathrm{\Gamma }_\alpha ,t_C)`$. In the following we discuss the physical picture of this approximation. The Green’s function $`G_{LL\sigma ,mn}^r`$ is a progapator (of the Hamiltonian $`H`$) for an electron from the left dot to the left dot. The first term $`g_{LL\sigma ,mn}^r`$ describes that the electron always stays in the left dot. The second term $`\underset{k_1}{}g_{LL\sigma ,mk_1}^rA_{R\sigma ,k_1n}g_{LL\sigma ,nn}^r`$ describes the progapation process in which the electron from left dot transits to the right dot while absorbing $`k_1k`$ photons, it then stays in the right dot for some time, then goes back to the left dot emitting $`nk`$ photons. The higher order term describes multi-time transit processes between the left dot and the right dot. In our approximation (e.g. for the 2nd term), we take $`k_1=n`$ and neglect others. We emphasize that this term is indeed the most important. For example, if $`k_1n`$, we consider an initial electron in the left dot sitting at the resonance state (e.g. at $`ϵ_L`$). This electron tunnels to the right dot absorbing $`k_1k`$ photons and goes back to the left dot emitting $`nk`$ photons. Then this electron has energy $`ϵ_L+(k_1n)\mathrm{}\omega `$, which has a difference $`(k_1n)\mathrm{}\omega `$ from the resonance state $`ϵ_L`$, and there is no state at energy $`ϵ_L+(k_1n)\mathrm{}\omega `$. Therefore this processes ($`k_1n`$) has very small probability to occur. In fact, when $`ϵ`$ is near a resonant state (e.g. $`ϵ_L`$), $`g_{\alpha \alpha \sigma }^r(ϵ)`$ \[Eq.(3) in the manuscript\] can be rewritten as $`1/(ϵϵ_L+i\delta )`$. Then the term in the left side of Eq.(C6) can be rewritten as:

$$\frac{A_{R\sigma ,k_1n}}{(ϵϵ_L+k_1\mathrm{}\omega +i\delta )(ϵϵ_L+n\mathrm{}\omega +i\delta )}$$
$`(C7)`$
where $`\delta `$ is a small positive real number ($`\mathrm{\Gamma }`$). Now it clearly shows that when $`k_1n`$, this term is always very small for any electron energy $`ϵ`$ \[note that $`A_{R\sigma ,k_1n}`$ has a small factor $`t_C^2`$\]. When $`k_1=n`$, this term is large at $`ϵ=ϵ_Ln\mathrm{}\omega `$. Therefore our approximation is reasonable if $`\mathrm{}\omega >>max(\mathrm{\Gamma }_\alpha ,t_C)`$.

Under this approximation, we can easily solve $`G_{LL\sigma ,mn}^r`$ from Eq.(C5)

$$G_{LL\sigma ,mn}^r(ϵ)=\frac{\delta _{mn}}{\left(g_{LL\sigma ,mm}^r\right)^1A_{R\sigma ,mm}}$$
$`(C8)`$
Then, $`G_{LR\sigma ,mn}^r`$ can be obtained immediately:

$$G_{LR\sigma ,mn}^r(ϵ)=G_{LL\sigma ,mm}^r(ϵ)\mathrm{\Sigma }_{LR\sigma ,mn}^rg_{RR\sigma ,nn}^r(ϵ)$$
$`(C9)`$
$`G_{RR\sigma ,mn}^r(ϵ)`$ and $`G_{RL\sigma ,mn}^r(ϵ)`$ can also be solved in the same manner.


format: text

text:
# Decomposing the 𝑆𝑈(𝑁) connection and the Wu-Yang potential

## I Introduction

Recently Faddeev et al. proposed a decomposition of the four dimensional $`SU(N)`$ Yang-Mills field $`A_\mu ^a`$ . In their paper with some ansatz the new variables were given for studying the knot theory and the QCD. One of our author (Duan) pointed out that the gauge potential should be decomposed in terms of the gauge covariant , i.e. $`A_\mu =a_\mu +b_\mu `$, which the $`a_\mu `$ satisfies the gauge transformation $`a_\mu ^{}=ga_\mu g^1+_\mu gg^1`$ and the $`b_\mu `$ satisfies the adjoint transformation $`b_\mu ^{}=gb_\mu g^1.`$ The $`a_\mu `$ part may show the geometry property of system and the $`b_\mu `$ part may be looked upon as vector boson which would be massive. Based on the decomposition theory of the $`SU(2)`$ gauge field , the gauge potential of the $`SU(N)`$ gauge field is decomposed in terms of local bases corresponding to the Cartan subalgebra without any hypothesis. With this decomposition the extended ’t Hooft electromagnetic tensor is derived and the Wu-Yang potential is given. At last, we discuss the sourceless solution of the gauge field equation of group $`SU(N).`$

## II SU(2) gauge field and the ’t Hooft monopole

In this section we introduce the decomposition theory of SU(2) gauge field and corresponding Wu-Yang potential. In terms of the $`SU(2)`$ gauge theory, the covariant derivation of a unit gauge field $`n^a(x)`$ is

$`D_\mu n^a(x)=_\mu n^a(x)+\epsilon ^{abc}B_\mu ^bn^c,a,b,c=1,2,3,`$

where $`B_\mu ^a`$ is the gauge potential of the $`SU(2)`$ gauge theory. By virtue of this definition we can give the expression of the gauge potential $`B_\mu ^a`$

$$B_\mu ^a=A_\mu n^a+\epsilon ^{abc}_\mu n^bn^c\epsilon ^{abc}D_\mu n^bn^c,$$
(1)
where $`A_\mu =B_\mu ^an^a`$ is the Abelian projection of the $`SU(2)`$ gauge potential. Thus the gauge potential is decomposed formally, one will find that the decomposition is very useful in the ’t Hooft monopole theory. We calculate the corresponding field strength

$`F_{\mu \nu }=_\mu B_\nu _\nu B_\mu [B_\mu ,B_\nu ],`$

from the Eq. (1) we give

$$F_{\mu \nu }^an^a=(_\mu A_\nu _\nu A_\mu )\epsilon ^{abc}n^a_\mu n^b_\nu n^c+\epsilon ^{abc}n^aD_\mu n^bD_\nu n^c.$$
(2)
Since the left hand $`F_{\mu \nu }^an^a`$ and the last term of right hand are the gauge invariant term, we find that the following quantity is a gauge invariant

$`f_{\mu \nu }=(_\mu A_\nu _\nu A_\mu )K_{\mu \nu },`$

where $`K_{\mu \nu }=\epsilon ^{abc}n^a_\mu n^b_\nu n^c`$. One can find that $`f_{\mu \nu }`$ is just the gauge field tensor defined by ’t Hooft which is fundamental for the magnetic monopole

$$f_{\mu \nu }=F_{\mu \nu }^an^a\epsilon ^{abc}n^aD_\mu n^bD_\nu n^c.$$
(3)
Let $`\stackrel{}{e}_a(x)`$ $`(a=1,2)`$ are the vierbein perpendicular to $`\stackrel{}{n}(x),`$ it is easy to verify that there exists a $`U(1)`$ potential

$$a_\mu =\epsilon ^{ab}(\stackrel{}{e}_a_\mu \stackrel{}{e}_b),$$
(4)
which satisfies $`K_{\mu \nu }=_\mu a_\nu _\nu a_\mu .`$ One can find that this potential is just the Wu-Yang potential of the magnetic monopole system.

## III Local basis and the SU(N) connection

Let $`T_a`$ ($`a=1,2,\mathrm{},r)`$ be Lie algebraic bases of the $`SU(N)`$ group $`G`$, and $`H_i(i=1,2,\mathrm{},L)`$ the Cartan subalgebra, i.e.

$$[T_a,T_b]=if_{abc}T_c,[H_i,H_j]=0.$$
(5)
The local basis of Cartan subalgebra is defined as

$$n_i(x)=U(x)H_iU^{}(x),$$
(6)
where $`U(x)`$ is a unitary matrix on manifold $`M.`$ The covariant derivative of the local basis $`n_i(x)`$ is

$$D_\mu n_i=_\mu n_iig[A_\mu ,n_i]$$
(7)
where $`A_\mu `$ is a $`su(N)`$ Lie algebra vector

$`A_\mu =A_\mu ^aT_a.`$

The curvature is defined as

$$F_{\mu \nu }=_\mu A_\nu _\nu A_\mu ig[A_\mu ,A_\nu ].$$
(8)
In terms of the relation

$$f_{abl}f_{acm}n_i^bn_j^c+n_i^ln_j^m=\delta ^{lm},$$
(9)
one can find that a $`su(N)`$ Lie algebra vector $`V`$ can be decomposed with the local basis

$$V=(V,n_i)n_i+[[V,n_i],n_i],Vsu(N).$$
(10)
Since it is a $`su(N)`$ Lie algebra vector, with the Eq. (7) the $`SU(N)`$ connect $`A_\mu `$ can be decomposed as

$$A_\mu =(A_\mu ,n_i)n_i+\frac{1}{ig}[_\mu n_i,n_i]\frac{1}{ig}[D_\mu n_i,n_i].$$
(11)
Similiar to the ’t Hooft electromagnetic tensor in the $`SU(2)`$ gauge field, we can define the extended ’t Hooft electromagnetic tensor

$$f_{\mu \nu }^i=(F_{\mu \nu },n_i)+\frac{i}{g}(n_i,[D_\mu n_j,D_\nu n_j]),$$
(12)
which is gauge invariant. By make using of the relation

$$(n_i,[n_j,A])=0,$$
(13)
one can prove that

$$f_{\mu \nu }^i=_\mu A_\nu ^i_\nu A_\mu ^i+\frac{i}{g}(n_i,[_\mu n_j,_\nu n_j]),$$
(14)
in which $`A_\mu ^i`$ is the projection of $`A_\mu `$ on the $`n_i`$

$`A_\mu ^i=(A_\mu ,n_i).`$

From the above discussion, one can find that this decomposition is different from the work of Faddeev et al. With the idea that the connection should possess its inner structure thus can be decomposed, we express the connection $`A_\mu `$ with the covariant part $`b_\mu `$ and the non-covariant part $`a_\mu `$, i.e. $`A_\mu =a_\mu +b_\mu `$. In the case of SU(N), we have

$$a_\mu =(A_\mu ,n_i)+\frac{1}{ig}[_\mu n_i,n_i],b_\mu =\frac{1}{ig}[D_\mu n_i,n_i].$$
(15)
In the section I, we discussed the non-covariant part $`a_\mu `$ should give the geometry information of system. In the following section we can find that this part corresponds to the Wu-Yang potential. In the paper of Faddeev et al., the decomposition of SU(N) was given by introducing the new variables and they focused on the covariant part $`b_\mu `$ to consider the effective Lagrangian.

## IV Wu-Yang potential and the sourceless solution of the SU(N) gauge field

Since it was proposed, the existence of Wu-Yang potential has been doubted for a long time. As a class of important gauge field configurations for QCD, Wu-Yang monopoles are fundamental for confinement and compete with the instanton-like configurations which are responsible for chiral symmetry breaking . In this section we will be sure that the Wu-Yang potential exist in the $`SU(N)`$ gauge theory.

From Eq. (6) one have

$$_\mu n_i=ig[\stackrel{~}{A}_\mu ,n_i],$$
(16)
where $`\stackrel{~}{A}_\mu =\frac{1}{ig}_\mu UU^{}`$, is obviously a flat connection, i.e.,

$$_\mu \stackrel{~}{A}_\nu _\nu \stackrel{~}{A}_\mu ig[\stackrel{~}{A}_\mu ,\stackrel{~}{A}_\nu ]=0$$
(17)
In terms of Eq. (14) and (17) we can obtain

$`{\displaystyle \frac{1}{ig}}(n_i,[_\mu n_j,_\nu n_j])=_\mu a_\nu ^i_\nu a_\mu ^i,`$

where

$$a_\mu ^i=(\stackrel{~}{A}_\mu ,n_i)$$
(18)
are the Wu-Yang potential. Thus we show that Wu-Yang potentials are Abelian projection of the flat connection in the $`SU(N)`$ gauge theory. Furthermore, since there are $`(N1)`$ Abelian projection directions, one can find $`(N1)`$ Wu-Yang potentials in $`SU(N)`$ gauge theory.

Finally, from Eq. (12) we note that when $`D_\mu n_i=0,`$

$$f_{\mu \nu }^i=(F_{\mu \nu },n^i).$$
(19)
But from the relation that

$$(D_\mu D_\nu D_\nu D_\mu )n_i=[F_{\mu \nu },n_i],$$
(20)
we find that $`[F_{\mu \nu },n_i]=0(i=1,2,\mathrm{},L).`$ That is to say, when $`D_\mu n_i=0,F_{\mu \nu }`$ must commute with $`n_i.`$ Therefore using Eq. (19), $`F_{\mu \nu }`$ can be expressed as

$$F_{\mu \nu }=f_{\mu \nu }^in_i.$$
(21)
Then from $`D_\mu n_i=0`$, we have

$`D_\nu F_{\mu \nu }=(_\nu f_{\mu \nu }^i)n_i.`$

Hence the solutions of equation

$$_\nu f_{\mu \nu }^i=0,i=1,2,\mathrm{},L,$$
(22)
corresponding to the solutions of the sourceless equation

$`D_\nu F_{\mu \nu }=0.`$

If we adopt the Lorentz condition $`_\nu A_\nu ^i=0(i=1,2,\mathrm{},L),`$ from Eqs. (14) and (22) one can get

$$^2A_\mu ^i=_\nu K_{\mu \nu }^i,$$
(23)
where

$`K_{\mu \nu }^i={\displaystyle \frac{i}{g}}(n_i,[_\mu n_j,_\nu n_j]).`$

The solution of Eq. (23) is

$$A_\mu ^i=\frac{1}{4\pi }G(xx^{})_\nu K_{\mu \nu }^i(x^{})d^4x^{},$$
(24)
where $`G(xx^{})`$ is the retarded Green’s function

$`^2G(xx^{})=4\pi \delta ^4(xx^{}).`$

Then, it follows from Eq. (11) that the equation

$$A_\mu =\frac{1}{4\pi }[G(xx^{})_\nu K_{\mu \nu }^i(x^{})d^4x^{}]n_i+\frac{1}{ig}[_\mu n_i,n_i],$$
(25)
is the sourceless solution of the gauge field equation.

## V Conclusion

In this paper, we give the decomposition of SU(2) and SU(N) connection with the idea that connection should possess their inner structures. One can find that the decomposition in this paper have no any hypothesis or ansatz since the expressions of connections are obtained with their definitions directly. In terms of the decomposition of SU(2) connection, we give the ’t Hooft invariant gauge field tensor which can describe the magnetic monopole correctly. And with the expression of SU(N) connection we consider the Wu-Yang potential which is responsible for the confinement in QCD.

## VI Acknowledgement

This work was supported by the National Natural Science Foundation of China and the Doctoral Foundation of the People’s Republic of China.


text:
# CO (J=7→6) Observations of NGC 253: Cosmic Ray Heated Warm Molecular Gas

## 1 Introduction

NGC 253 is a nearby ($`d2.5`$ Mpc; Mauersberger et al., 1996 ), nearly edge-on ($`i78^{}`$, Pence 1981), isolated starburst galaxy with the starburst strongly concentrated at the nucleus. Based on the distribution of the 10–30 $`\mu `$m continuum, the bulk of NGC 253’s starburst is confined to a $`60`$ pc region centered just south-west of the dynamical nucleus (Telesco, Dressel, & Wolstencroft, 1993). The 1–300 $`\mu `$m luminosity within about 30′′ is $`1.6\times 10^{10}\mathrm{L}_{}`$ (Telesco & Harper, 1980) while the total IR luminosity of NGC 253 detected by IRAS is $`\mathrm{L}_{\mathrm{IR}}2\times 10^{10}\mathrm{L}_{}`$ (Rice et al., 1988). There is a clear bar in the 2MASS image of the galaxy, and the kinematics show evidence for orbits in a bar potential (Scoville et al. 1985; Das et al. 2001, and references therein). The bar likely plays a role in transporting gas to the nucleus and thus supporting the starburst. Current estimates of the neutral gas mass in the central $`20^{\prime \prime }50^{\prime \prime }`$ range from $`2.5\times 10^7\mathrm{M}_{}`$ (Harrison, Henkel, & Russell, 1999, hereafter HHR99) to $`4.8\times 10^8\mathrm{M}_{}`$ (Houghton et al., 1997).

As is the case in many of the nearby starburst galaxies, much of the molecular gas in NGC 253 appears to be highly excited (Wild et al., 1992; Mao et al., 2000; Ward et al., 2002). Observations of $`\mathrm{J}=43`$ and $`\mathrm{J}=65`$ transitions of CO (Israel & Baas, 2002; Harris et al., 1991), as well as $`\mathrm{J}=10`$ and $`\mathrm{J}=32`$ transitions of HCN (Paglione et al., 1995, 1997) suggest that a large fraction of the molecular gas in the nucleus of NGC 253 has $`\mathrm{n}_{\mathrm{H}_2}>10^4\mathrm{cm}^3`$ and $`\mathrm{T}100\mathrm{K}`$ or higher. Observations of higher-J CO transitions are critical to measure this warm molecular gas component, and probe its physical conditions. In our Galaxy, the warm, dense molecular gas conditions which excite the mid-J CO lines (J = 4–8) are typically associated with molecular cloud surfaces exposed to stellar far-UV radiation. In these photodissociation regions (PDRs), the UV dissociates and ionizes the gas as it enters the cloud and is absorbed by dust. The result is a layered geometry based on the extinction of UV; the progression from the exterior inward is ionized, neutral atomic, warm molecular, and finally cool molecular gas on the interior (Tielens & Hollenbach, 1985; Wolfire, Tielens, & Hollenbach, 1990; Kaufman et al., 1999). Several of these regions have been studied in star formation sites within the Galaxy (Harris et al., 1985; Schmid-Burke et al., 1989; Stacey et al., 1993), and they have been discussed in the context of starburst nuclei (Mao et al., 2000; Carral et al., 1994). Other potential molecular gas heating sources are outflows and cloud-cloud collisions in which shocks convert dynamical energy to thermal energy of the gas, such as in the circumnuclear disk of the Galactic Center, and in the massive star formation sites W 51 and DR 21 (Jaffe, Harris, & Genzel, 1987; Jaffe et al., 1989). Finally, it has been proposed that much of the molecular gas heating in starburst nuclei may be due to cosmic rays (Suchkov, Allen, & Heckman, 1993, hereafter SAH). To better understand the energy balance of the molecular gas in the starburst nucleus of NGC 253, we have observed the CO $`\mathrm{J}=76`$ transition with our submillimeter imaging spectrometer, SPIFI on the JCMT.

## 2 Observations

NGC 253 was observed in April 1999 during the commissioning run of the South Pole imaging Fabry-Perot interferometer (SPIFI) (Bradford et al., 2002) at the James Clerk Maxwell Telescope (JCMT) on Mauna Kea, Hawaii. The Fabry-Perot interferometer (FPI) was set to a velocity resolution of 70 km s-1 and was scanned over 620 km s-1 with the center wavelength of the scan set to the redshifted (250 km s-1) CO $`\mathrm{J}=76`$ rotational transition. The FPI scan length was measured with interference fringes of a He-Ne laser, and absolute frequency reference was an absorption spectrum of the $`\mathrm{J}=67`$ CO transition from gas in a cell against a hot load. On this observing campaign the bolometer array was equipped with 12 functioning detectors arranged on a $`4\times 4`$ grid with four blanks. Because of a figuring error in the polyethylene lenses used as a f-converter, the main beam as measured with Mars had a size of $`11.5^{\prime \prime }`$ (FWHM), 1.6 times bigger than the diffraction-limited beam size of $`7^{\prime \prime }`$ per pixel. Mars was also used for flux calibration assuming a brightness temperature of 201 K (Wright, 1976; Rudy et al., 1987). Flat-fielding of the array was accomplished by rastering Mars around the array. For the period of the observations, the zenith atmospheric transmission was between 0.22 and 0.29, giving a transmission to the source between 0.05–0.12. This low and variable atmospheric transmission provided the greatest calibration uncertainty, which we estimate to be about 30%. For the observation of NGC 253 the array was centered at the IR nucleus at $`R.A.=0^h45^m06^s.0`$, $`Dec.=25^{}33^{}40^{\prime \prime }`$ (1950) ((0,0)-position). We estimate our pointing accuracy to be 5′′. The data were obtained by chopping and nodding the telescope and the data presented here were obtained with 20 minutes on-source integration. The 50 minute duration of the observations resulted in a field rotation of 18 about the center of the array, producing a total travel of 7′′ ($`\pm `$3.5′′ about the mean) for the corners of the array. Because this is comparable to the pointing uncertainty, and is substantial only for the corner pixels, we neglect the rotation and adopt the mean position angle.

## 3 Results

### 3.1 Comparison with Other Molecular Tracers

Figure 1 shows the spectra obtained in this single array pointing toward the nucleus of NGC 253. The data are plotted in main beam brightness temperature units, correcting for beam efficiency, atmospheric transmission, and source coupling. Typical spectra peak at 3–5 K, with a single spectrum showing a peak of 12 K; integrated intensities are between 400 and 1200 K km s-1. The integrated intensities in K km s-1 are roughly half of the CO $`\mathrm{J}=65`$ intensities of Harris et al. (1991), when the $`\mathrm{J}=65`$ antenna temperatures are corrected to a $`10^{\prime \prime }`$ main beam. The peak temperature of 12 K at $`v=150`$ km s-1 in the spectrum just northeast of the nucleus has a similar peak above the Gaussian profile in the $`\mathrm{J}=65`$ spectrum toward the nucleus. In addition to the bright spectra near the nucleus, we detect emission at 30–40% of the nuclear value at positions 10′′–20′′ offset from the nucleus. Harris et al. also detect strong CO $`\mathrm{J}=65`$ emission at offset positions. They measure an integrated intensity 60% of the nuclear value as far as $`14^{\prime \prime }`$ offset, consistent with our observations.

Interferometric observations of millimeter wave CS (Peng et al., 1996), CO (Canzian, Mundy, & Scoville, 1988; Carlstrom et al., 1990), and SiO (Garcia-Burillo et al., 2000) transitions show evidence for a nuclear bar with major axis oriented 70 east of north. The velocity field of this bar shows a 100$`\mathrm{km}\mathrm{s}^1`$ shift from the northeast to the southwest. Though the $`\mathrm{J}=76`$ map does not have the resolution or sampling to spatially resolve the bar, it also shows the velocity shifting from 150-200 $`\mathrm{km}\mathrm{s}^1`$ in the northeast to 250-300 $`\mathrm{km}\mathrm{s}^1`$ in the southwest, showing the overall kinematics of gas in the bar, in good agreement with the interferometric observations. The observed line widths are between 200 and 300$`\mathrm{km}\mathrm{s}^1`$, consistent with the 150 to 240 $`\mathrm{km}\mathrm{s}^1`$ line widths observed in the lower-J CO lines with 12 to 15′′ beams (Mauersberger et al. (1996); HHR99), convolved with our 70 $`\mathrm{km}\mathrm{s}^1`$ Lorentzian instrumental profile.

### 3.2 Model fit to the nuclear excitation conditions

To make quantitative estimates of the physical parameters of the nuclear molecular gas we have compiled all CO rotational transitions observed toward the central 15′′ available in the literature to compare with predictions of an escape probability, large velocity gradient model for CO line excitation. In contrast with previous studies of the molecular material in starburst galaxies which invoke two components (HHR99) or complex PDR models (Mao et al., 2000), our modeling suggests that the observed $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ intensities can be reproduced with the bulk of the gas mass at fairly high excitation. The line ratios can be fit with a range in conditions from $`\mathrm{T}=100\mathrm{K}`$, $`\mathrm{n}_{\mathrm{H}_2}5\times 10^4\mathrm{cm}^3`$ to $`\mathrm{T}>300\mathrm{K}`$, $`\mathrm{n}_{\mathrm{H}_2}2.0\times 10^4\mathrm{cm}^3`$. Recent measurement of rotational H2 lines (Rigopoulou et al., 2002) suggests that the temperatures may be on the low end of the allowed range, and we thus we favor $`\mathrm{T}100\mathrm{K}`$ and $`\mathrm{n}5\times 10^4\mathrm{cm}^3`$. Combining the 13CO transitions constrains the column density to be only slightly larger than that measured by HHR99 in a 23′′ beam. To reproduce the observed intensities with the LVG formalism requires large velocity gradients ($`\mathrm{dv}/\mathrm{dz}1080\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$), which are much higher than is traditionally used, but nevertheless plausible for the starburst nucleus. Here we describe our modeling approach and the basis for these conclusions.

Including our results, 12CO has been observed toward the nucleus of NCG 253 in all transitions up to $`\mathrm{J}=76`$ except for $`\mathrm{J}=54`$, with a variety of beam sizes. Also, the optically thin isotopes 13CO and C18O have been observed in transitions up to $`\mathrm{J}=32`$. Table 1 summarizes the 12CO and 13CO observations and the corrections to refer them to the 15′′ beam. Measurements for all $`\mathrm{J}=32`$ and lower transitions are from HHR99. Coupling corrections, where necessary, are based on power laws. The $`\mathrm{J}=21`$ lines have been observed with both 12′′ and 23′′ beams, and the resulting power law is $`\mathrm{I}\theta ^{0.76}`$. For $`\mathrm{J}=10`$, there are no single dish observations with a beam smaller than 23′′, but based on the centrally peaked $`\mathrm{J}=32`$ / $`\mathrm{J}=21`$ ratio map, the $`\mathrm{J}=10`$ emission is probably more extended than $`\mathrm{J}=21`$. We adopt a power law of $`\mathrm{I}\theta ^{0.4}`$, but estimate a 20% uncertainty in the derived intensity, which is sufficient to accommodate a power law exponent between 0 and 0.76. For $`\mathrm{J}=32`$, the intensities are measured from the published spectra in 15′′ beams. For $`\mathrm{J}=43`$ there are published measurements with both a 10.6′′ beam (Israel et al., 1995) and a 22′′ beam (Israel & Baas, 2002). Allowing these two measurements to define a power law gives $`\mathrm{I}\theta ^{1.02}`$, and we use this to determine the intensity within 15′′. The $`\mathrm{J}=65`$ measurement is from Harris et al. (1991) in $`\mathrm{T}_\mathrm{R}^{}`$ units, and is corrected to 15′′ according to their stated 8′′/30′′ composite beam profile. The $`\mathrm{J}=76`$ observation from this work is referred to an 11.5′′ main beam based on observations of Mars, with a small correction for the 15.8′′ planetary disk. The estimated uncertainties are also included in Table 1. The beam filling corrections become the dominant source of uncertainty in the $`\mathrm{J}=10`$ (described above) and $`\mathrm{J}=21`$ lines, we estimate this to be about 15%. For the $`\mathrm{J}=32`$ transition in a 15′′ beam, we estimate 20% for the overall uncertainty including systematics and statistics. For the $`\mathrm{J}=43`$ and higher submillimeter transitions, the uncertainties in the measurements themselves are substantial, as much as 30%.

The LVG model calculates intensities of rotational CO lines given a molecular hydrogen density, a temperature, and a velocity gradient. The calculation is similar to that used in other multi-line analyses, iteratively solving for the level populations up to $`\mathrm{J}=15`$. We use the escape probability $`\beta =\left(1e^\tau \right)/\tau `$ of Goldreich, & Kwan (1974) appropriate for spherical clouds undergoing gravitational collapse. H2 – CO collisional excitation rates are those of Flower (2001), increased by 21% to account for collisions with He (McKee et al., 1982; Viscuso & Chernoff, 1988; Flower & Launay, 1985; Flower, 2001; Schinke et al., 1985).

#### 3.2.1 12CO Intensity Ratios

A sensible approach for comparing an LVG model with observations is to use the line intensity ratios to determine the physical parameters (temperature, density, and opacity) while the absolute line intensities determine the total column density. In addition to the 12CO intensity ratios, the optically-thin isotopes can be useful in constraining the 12CO optical depths, which are greater than unity for the lines of interest. In order to concentrate on the nuclear gas, we do not put much emphasis on the $`\mathrm{J}=10`$ transition. Since it connects the ground state and and the $`\mathrm{J}=1`$ level is only 6 K above ground, the $`\mathrm{J}=10`$ transition traces any molecular material in the beam. A related and more confusing problem is the potential for self-absorption by quiescent material along the line of sight but exterior to the nucleus, a particular concern for nearly edge-on systems such as NGC 253. This may be the reason that any model which fits the $`\mathrm{J}=21`$ and higher transitions over predict $`\mathrm{J}=10`$. In an approach similar to that of HHR99, we invoke a small diffuse component which is negligible energetically, but which emits primarily in $`\mathrm{J}=10`$ to account for its observed intensity.

For a given LVG calculation, the parameters which determine the line intensity ratios are the temperature, density, and column density of CO per unit velocity interval, $`\mathrm{N}_{\mathrm{CO}}/\mathrm{dv}`$, called the opacity parameter. $`\mathrm{N}_{\mathrm{CO}}/\mathrm{dv}`$ is a combination of density, velocity gradient of the gas, and CO abundance, $`\mathrm{N}_{\mathrm{CO}}/\mathrm{dv}=\mathrm{X}_{\mathrm{CO}}\mathrm{n}_{\mathrm{H}_2}(\mathrm{dv}/\mathrm{dz})^1`$ where $`\mathrm{dv}/\mathrm{dz}`$ is the velocity gradient in $`\mathrm{km}/\mathrm{s}\mathrm{pc}^1`$. We have assumed a CO abundance $`\mathrm{X}_{\mathrm{CO}}`$ of $`8\times 10^5`$, the value measured for galactic clouds (Freking, Langer, & Wilson, 1982). While the metallicity in the nucleus of NGC 253 is as much as 3 $`\times `$ solar (Mauersberger et al., 1996), the large observed \[CI\] intensities suggest that, relative to the Galaxy, the gas-phase carbon budget includes a much larger neutral carbon fraction, and smaller CO fraction (Israel & Baas, 2002), thus the CO abundance is not clearly much higher in spite of the large metallicity. Intensity ratio calculations for a higher CO abundance are identical to those with the canonical $`\mathrm{X}_{\mathrm{CO}}`$ value, but with the velocity gradient decreased by the same ratio. Raising the CO abundance or lowering the velocity gradient increases the the optical depth in the lines. Figure 2 plots the results of LVG calculations of diagnostic line ratios as a function of gas temperature, density, and velocity gradient. We have found that the $`\mathrm{J}=76`$ to $`\mathrm{J}=32`$ intensity ratio is a good measure of the overall excitation, and that matching this ratio typically satisfies all the lines above $`\mathrm{J}=10`$ to within the observational uncertainties.

As can be seen from the solid contours in Figure 2, the 12CO $`\mathrm{J}=76`$ to $`\mathrm{J}=32`$ intensity ratio of 7.7 means that the gas is of high excitation. Because the run of line intensity with J is increasing up to $`\mathrm{J}=76`$, there is no unique solution, nor even upper limits to the temperature or density. Lower limits, however, are strong. For temperatures less than 300 K, the density must be greater than $`23\times 10^4\mathrm{cm}^3`$. A temperature less than 150 K requires a density greater than $`4\times 10^4\mathrm{cm}^3`$, clearly different from the typical Galactic molecular cloud. To first order, these results are independent of the velocity gradient and CO abundance because the lines are optically thick, even at velocity gradients as large as 80 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$.

#### 3.2.2 $`{}_{}{}^{13}\mathrm{CO}`$ Intensities

If their abundances relative to the main isotope are known, the intensities in the 13CO transitions can be used to determine the 12CO optical depths. This is the approach taken by HHR99 in their study of the nuclear regions of NGC 253 with the lower-J lines. The key to their analysis is the assumption of a 12CO to 13CO abundance ratio of 40–50, originating from analysis by Henkel et al. (1993); Henkel & Mauersberger (1993). Our results are consistent with a ratio on the upper end of this range, and we have chosen 50 as a value to use in the model plots. Included in Figure 2 as dashed lines are the 12CO to 13CO $`\mathrm{J}=32`$ intensity ratios predicted by the LVG models. The observed ratio of 12 is plotted as a thick dashed line. To correctly obtain the line intensity ratios between the two isotopes requires running the full LVG model for the optically thin species, since radiative trapping is important in determining the level populations of 12CO but is much less so for the optically thin species. Because the 13CO transitions are optically thin for all relevant conditions, the values plotted in Figure 2 can be shifted to accommodate a different isotope ratio. In the density-temperature plane, a larger isotope ratio is considered by choosing a line ratio contour with a value decreased from the observed 12 by the same factor.

In the panels of Figure 2, the temperature and density at which the two observed line ratios meet is the regime of possible solutions. It is evident that small velocity gradients are unable to generate the observed line ratios. At these small velocity gradients, the conditions which excite upper-J lines generate large optical depths in the 12CO lines, and their fluxes are too small compared with the optically thin isotope. Increasing the velocity gradient does provide possible solutions, again, with high excitation conditions in both temperature (T$`>`$100 K) and density ($`\mathrm{n}_{\mathrm{H}_2}>4\times 10^5`$). Though unconventional in LVG analyses, the large velocity gradient is plausible for the starburst nucleus, and is discussed in Section 3.2.4.

#### 3.2.3 CO Column Density

Even though the 12CO lines are optically thick, their observed intensities put lower limits on the CO column density independent of the optically thin isotopes. The intensity in a line is determined by the product of column density in the upper level and the optical depth factor $`1e^\tau `$, hence a lower optical depth (larger velocity gradient) generates more emission for a given CO column density. The curves in Figure 3 show the column density required by the LVG model to reproduce the observed intensities. For each velocity gradient, the curve is parametrized in terms of temperature, with the density determined according to the observed 12CO line ratio shown in Figure 2. Examination of Figure 3 indicates that the column density of warm gas must be substantial. For example, if the gas has a temperature below 200 K and velocity gradient below 40 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, the column density must be larger than $`3\times 10^{18}\mathrm{cm}^2`$, a lower limit which is comparable to the $`3.5\pm 1.5\times 10^{18}\mathrm{cm}^2`$ estimated in the central 23′′ with the $`\mathrm{J}=32`$ and lower transitions (HHR99).

To improve upon the lower limits provided by the 12CO transitions alone, we incorporate the 13CO intensities. The boxes in Figure 3 show the column densities required by the LVG model under the conditions corresponding to the intersections of the two diagnostic line ratio curves in the panels of Figure 2. Boxes are plotted for values of $`\mathrm{dv}/\mathrm{dz}`$ = 10, 20, 40 and 80 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, as velocity gradients less than 10 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$ do not allow reasonable solutions. The sizes of the boxes correspond to a $`\pm `$30% uncertainty in temperature, and an corresponding uncertainty in column density given by local slope of the T-$`\mathrm{N}_{\mathrm{CO}}`$ relation. The locus of these regions indicates that the CO column density is between 3–7$`\times 10^{18}\mathrm{cm}^2`$, for all allowed values of the excitation. This is consistent with the 12CO upper limits, and is comparable to the value derived by HHR99 for the central 23′′.

The results of our LVG analysis can be summarized as follows: the molecular gas in the central 180 pc of NGC 253 is much more highly excited than that traced with only the lower-J lines in a larger 280 pc beam. The gas is warmer than 100 K, and has a density greater than $`2\times 10^4\mathrm{cm}^3`$. The 12CO lines have optical depths not larger than 2–4, requiring velocity gradients $`>10`$ $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, larger than traditionally employed in LVG studies. Even larger velocity gradients ($`>20`$ $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$) are required if the temperature is constrained to be below 200 K. The CO column density is well-determined by the 13CO lines: $`\mathrm{N}_{\mathrm{CO}}=`$3–7$`\times 10^{18}\mathrm{cm}^2`$, only slightly larger than the values estimated in a larger beam. Finally, while the 13CO $`\mathrm{J}=10`$ transition requires an additional component of cool and/or diffuse gas, this component is small by mass, negligible energetically, and may not even be associated with the nuclear region. There is no reason to invoke a large component of cool gas as would be required, for example, in a study of the central regions of the Milky Way.

#### 3.2.4 The Large Velocity Gradient

Previous LVG studies of NGC 253 and M82 have employed a velocity gradient of around 1 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, which is approximately the ratio of the total observed velocity width to the size of the observed region (HHR99; Mao et al. (2000)), and comparable to the values used in studies of Galactic molecular clouds. This is a lower limit to the true velocity gradient in the material, appropriate if it were only subject to the central potential of the galaxy with no local motion. In reality, clouds and clumps will be subject to their own gravity as well as that of the stars around them. The escape velocity of a self gravitating cloud of radius r, and density $`\rho `$ determines a velocity gradient given by $`\mathrm{dv}/\mathrm{dr}\sqrt{8\pi /3\mathrm{G}\rho }`$, which gives 10 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, for a density of $`\mathrm{n}_{\mathrm{H}_2}4.5\times 10^4\mathrm{cm}^3`$. Models on the low end of the allowed temperature range, e.g. T $``$ 120 K require a velocity gradient of order 40 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, a factor of 4 larger than this escape velocity. This is not implausible; if the star-forming clouds in NGC 253 are similar to the Orion molecular cloud core, then one might expect the stellar mass density to be about 4 times the gas density (O’Dell, 2001), which would double the effective escape velocity. Furthermore, Carral et al. (1994) find the pressure of the ionized gas in the central 45” of NGC 253 is $`7\times 10^6\mathrm{K}\mathrm{cm}^3`$ \- larger than the thermal pressure within our modeled molecular gas. This high pressure ionized gas has a large volume filling factor (0.01 to 0.3) within their much larger (45′′) beams. It is thus quite plausible that within our 15′′ modeled region the ionized gas place a important role in cloud confinement, further increasing the velocity dispersion in the gas. In short, the large velocity gradients required for a good model fit are very plausible given the dynamics in the nucleus of NGC 253.

### 3.3 Comparison with rotational H2 lines

Because the $`\mathrm{J}=76`$ transition is not high enough to see the turn-over in CO excitation, the CO analysis does not provide a good upper limit to the gas temperature. Fortunately, we can use the ISO H2 pure rotational line observations to provide an important temperature constraint. Rigopoulou et al. (2002) have recently reported ISO SWS observations of the pure rotational lines of H2 from the central 14′′ $`\times `$ 20′′ regions of NGC 253. These quadrupole H2 lines have low radiative rates, so they are easily thermalized and optically thin. In principle the lines can be used to measure temperature and column density. However, even the lowest H2 rotational transition energy is 515 K, thus the observed line emission is heavily biased toward the very warmest regions (which may be negligible by mass), and each observed transition is an exponential function of gas temperature. Not surprisingly, the line intensities are not consistent with a single excitation temperature. Based on the S(0), S(1), S(2), S(5), and S(7) lines, Rigopoulou et al. (2002) conclude that there must be a substantial mass of warm, dense gas in the nucleus of NGC 253, broadly consistent with our results. More specifically, the modeled molecular hydrogen column density of $`\mathrm{N}_{\mathrm{H}_2}=\mathrm{N}_{\mathrm{CO}}/\mathrm{X}_{\mathrm{CO}}`$4–8$`\times 10^{22}\mathrm{cm}^2`$ is consistent with the S(0) data presented by Rigopoulou et al. (2002) for gas temperatures in the range of 95–120 K. We therefore favor models on the low end of the temperature range allowed by CO analysis. Because of the exponential dependence on temperature, this H2 temperature constraint is a rather weak function of H2 column density, hence also a weak function of CO column density and CO abundance.

We combine this H2 result with the parameter ranges allowed by the LVG analysis, and the physical constraint of a small velocity gradient to arrive at our preferred model, shown in Figure 4. This model uses $`\mathrm{dv}/\mathrm{dz}`$ = 40 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, toward the bottom of the range allowed by the low temperature and the observed 13CO intensity. As described above, the lower velocity gradients are more consistent with gravitationally bounds clouds. The modeled conditions are $`\mathrm{T}=120\mathrm{K}`$, $`\mathrm{n}_{\mathrm{H}_2}=4.5\times 10^5\mathrm{cm}^3`$, and the CO column density required to generate the observed line intensities is $`5.3\times 10^{18}\mathrm{cm}^2`$. At this temperature, about 40$`\%`$ of the observed H2 S(1) flux is generated in the bulk component, and the contributions to the higher-J H2 lines (which require excitation of levels with $`\mathrm{E}>1500\mathrm{K}`$) are less than 10% of the observed line fluxes. As an estimate of the uncertainties, we note that even with the temperature constrained, a larger velocity gradient still allows a smaller column density. Decreasing the velocity gradient, or increasing the CO abundance can permit a larger column density, but this quickly becomes inconsistent with the 13CO observations unless the 13CO relative abundance is decreased below our adopted 1/50. Combining all the uncertainties in the abundance ratio, modeling, and observations, we estimate an uncertainty of 40% for the final adoped column density: $`\mathrm{N}_{\mathrm{CO}}3.57\times 10^{18}\mathrm{cm}^2`$. Our model parameters are summarized in Table 2.

### 3.4 Warm Gas Mass

With the assumed CO abundance relative to H2 of $`8\times 10^5`$, the CO column density corresponds to a mass of excited molecular gas in the central 15′′ of 2.5–5.1$`\times 10^7\mathrm{M}_{}`$. In reality, the CO abundance is not well-known. The abundance may be larger than the canonical Galactic value, given the 3$`\times `$ solar metallicity (but see discussion above), and this would result in a smaller gas mass by as much as a factor of 3. Based on an assumed carbon gas to dust ratio and the measured fractions of gaseous carbon bearing species, HHR99 infer a CO abundance relative to H2 of $`2.2\times 10^4`$, that is, about the factor of 3 larger than our adopted value. With this abundance, the gas mass would be smaller: 0.9–1.8$`\times 10^7\mathrm{M}_{}`$. We do note a larger CO abundance demands even larger velocity gradients (more than 30 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$) to allow working solutions. For comparison, the most recent estimate of the gas mass based on 1.3 mm dust continuum is $`4.4\times 10^7\mathrm{M}_{}`$ (HHR99, Krügel et al. (1990)). Since the 1.3 mm continuum has a half-power size of only 16′′ $`\times `$ 11′′ much of this flux arises within our within our 15′′ beam, and traces the same material. The agreement with our derived gas mass suggests that while our CO abundance may be slightly high, it is not off by more than a factor of 2.

### 3.5 Suitability of Multiple-Component Models

It is traditionally thought that mid-J CO emission originates in the PDR surfaces of molecular clouds exposed to UV starlight. In this scenario, the warm excited gas is separate from the bulk of the molecular gas, and usually represents a small fraction of the total mass. Observations of several PDR regions in the Galaxy have motivated modeling which includes detailed PDR chemistry and radiative transfer to predict the line emission. Köster et al. (1994) have constructed detailed models of PDRs on the surfaces of small clumps to predict 12CO and 13CO line intensities in star-formation regions. Relative to a simple plane-parallel model, the clumpy geometry allows an increase in local gas density, and allows UV photons to penetrate to a greater fraction of the gas, increasing the overall molecular excitation level. While such models are successful in reproducing the intensity ratios between 12CO lines, they fail to account for both the large integrated intensities and the large 12CO to 13CO intensity ratios observed in a starburst nucleus such as NGC 253. For example, Köster et al. (1994) model slabs with thickness given by $`\mathrm{A}_\mathrm{V}=10`$, illuminated from both sides, including chemistry of both 12C- and 13C-bearing species. The ratios among the 12CO lines are well-fit with $`\mathrm{n}_{\mathrm{H}_2}=10^5\mathrm{cm}^3`$, $`\chi _{\mathrm{UV}}=10^5`$ Habing ($`=1.6\times 10^2\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$), but to generate the observed $`\mathrm{J}=76`$ flux in the central 15′′ of NGC 253 with a superposition of these regions would require between 5 and 20 times the CO column density estimated with our LVG analysis. The 13CO lines show the same effect – the predicted $`\mathrm{J}=32`$ 12CO to 13CO intensity ratio is around 3, compared with the observed value in 15′′ of $`12`$ (HHR99). Thus the existing PDR models are not appropriate for the global properties of a starburst nucleus. The key problem is the high optical depth in the PDR models, resulting from the small assumed velocity width per gas column. The effective velocity gradient of the clumpy PDR model is around 9 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$, too low as we have discovered in Section 3.2.

To consider a PDR scenario with the appropriate optical depths in the lines for NGC 253, we have examined LVG models with two distinct gas components. Though this approach is not as rigorous as a full treatment of the chemistry and UV penetration used by the Galactic PDR models, to first order it accounts for the combination of warm UV-exposed surfaces and cool cloud interiors which together would generate the observed CO emission. For NGC 253, we have found that while it is possible to fit the observed intensities without the full mass of gas being highly excited, the warm fraction must be substantial. To reproduce the observed $`\mathrm{J}=65`$ and $`\mathrm{J}=76`$ transitions, at least 30% of the total column density estimated with the 13CO ratio must have high excitation conditions. Furthermore, the curves in Figure 3 shows that this scenario requires an extremely large velocity gradient ($``$100 $`\mathrm{km}\mathrm{s}^1\mathrm{pc}^1`$) and/or very high temperatures ($``$300 K) in this 1/3 warm fraction in order to reproduce the observed intensities. Moreover, as we outline in Section 4.1.1, even with only 1/3 of the total column density accounting for the high-J line emission, the amount of warm molecular gas would still be 7–10$`\times `$ greater than the PDR-associated atomic gas, so the UV heating is not a viable origin of the observed CO emission.

## 4 Discussion

### 4.1 Energetics of the Molecular Gas in NGC 253

Independent of the details of the excitation and radiative transfer model, the general shape of the CO intensity vs. J curve measures the luminosity of the molecular gas in the central 15′′ of NCG 253. Taking our favored model plotted in Figure 4 as a guide, the total flux in the CO lines from the modeled region is $`1.5\times 10^{11}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2`$. This is 25% of the \[C II\] 158 $`\mu `$m and \[O I\] 63 $`\mu `$m intensities from the entire inner $`45^{\prime \prime }`$ (Carral et al., 1994), and a substantial $`1.4\times 10^4`$ fraction of the total IR luminosity of the galaxy. At a distance of 2.5 Mpc, this molecular gas flux corresponds to a luminosity of $`2.9\times 10^6\mathrm{L}_{}`$, giving a luminosity to mass ratio of $`0.1L_{}/M_{}`$ for the gas containing CO based on the mass estimates of Section 3.4. The large luminosity in the CO lines is a simple consequence of the fact that such a large fraction of the gas mass is highly excited. In what follows, we examine potential heating sources to generate the large luminosity. Simple arguments show that traditional heating sources for warm molecular gas – far-UV photons and cloud-cloud collisions – are unlikely to contribute substantially. In contrast, the enhanced cosmic ray flux of the starburst nucleus can easily match the observed cooling in the CO lines, and is likely responsible for the high excitation conditions inferred for the bulk of the molecular ISM.

For completeness we note that the CO lines are a good measure of the total cooling of the molecular gas. The five H2 lines observed by Rigopoulou et al. (2002) suggest that all the rotational lines combine to produce a flux of at most $`1\times 10^{11}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^2`$, less than the CO lines. Another potential coolant is C0, which can have an abundance comparable to that of CO. However, it only has two transitions each of which could be comparable to a single CO line if optically thick and thermalized; thus C0 cannot be a dominant coolant. Also, we note that though the modeling described above carries uncertainties as great as 40% in the derived parameters, the discussion below typically concerns factors of 3–30, and the general arguments hold independent of the modeling details.

#### 4.1.1 Stellar UV Heating

In Galactic starformation regions, warm CO is typically produced by UV impinging on the surfaces of molecular clouds. As it penetrates a cloud, the UV is attenuated, producing in turn ionized, atomic, and warm molecular layers before it is extinguished (Tielens & Hollenbach, 1985). The thermal pressure implied by our model fit to the CO lines is $`5.4\times 10^6\mathrm{K}\mathrm{cm}^3`$, similar to the pressure derived with atomic PDR line ratios of $`2.9\times 10^6\mathrm{K}\mathrm{cm}^3`$ (Carral et al., 1994). This is consistent with a model in which the PDR gas forms the surface of the warm molecular clouds. The penetration of UV into a cloud is determined by the dust opacity, proportional to the gas column density according to

$$A_\mathrm{V}=2N_{\mathrm{H}_2}/1.9\times 10^{21}\mathrm{cm}^2.$$
(1)
Here we have assumed solar dust abundance. (This is not unreasonable—though the metallicity in the nucleus of NGC 253 may be 3$`\times `$ solar (Mauersberger et al., 1996), there is evidence for substantial destruction of dust grains in an unusually large gas-phase Si abundance (Carral et al., 1994)). In standard Galactic PDRs, the first $``$ 3 magnitudes from the surface of the cloud are photodissociated and cool through the atomic lines, while the warm molecular gas exists from $`\mathrm{A}_\mathrm{V}`$ 3–6. To first order, the masses of these two components are comparable. Yet in the nucleus of NGC 253, Carral et al. (1994) measure $`2.4\times 10^6`$M of photodissociated gas traced by the 158$`\mu \mathrm{m}`$ \[CII\] line in a large (45′′) beam. Scaling by the observed molecular and far-IR continuum morphology (Canzian, Mundy, & Scoville, 1988), (Telesco, Dressel, & Wolstencroft, 1993), we estimate only 56% of the photodissociated gas traced with the large beam falls within our 15′′ beam. Therefore, the photodissociated gas mass is about a factor of 20–30 smaller than the 2–5$`\times 10^7`$ M of warm molecular gas. There is no reasonable geometry in which two layers of comparable thickness can have volumes different an order of magnitude. It is possible that the distribution of heating between the molecular and atomic components could be made to favor the molecular gas with soft or weak UV fields relative to standard PDR templates. However, the conditions derived with the atomic lines are consistent with standard photodissociation region (PDR) models – the UV field can be understood as arising from O 7.5 stars with an effective temperature of 34,500 K (Carral et al., 1994; Engelbracht et al., 1998), so scenarios invoking soft or weak UV fields are unlikely.

As pointed out in Section 3.5, as little as 30% of the total CO could provide the line luminosities in a warm component, emitting more efficiently lower opacity in the lines. The remaining 70% of the molecular gas would then be in a cool, more diffuse component, possibly the UV-shielded cloud cores. This scenario is attractive because UV-shielded regions are required to host the HCN observed toward the nucleus (Paglione et al., 1995), which has a low dissociation energy ($``$ 5 eV). However, a factor of 3 reduction in warm CO-containing gas mass would still leave a factor of 7–10 more warm molecular gas than atomic gas. Thus standard PDR models which invoke stellar far-UV photons to heat the surface of molecular clouds are not a good match for the heating of the bulk of the molecular gas in the nucleus of NGC 253.

#### 4.1.2 Dynamical Heating

Another potential explanation for the large molecular luminosity is dynamical energy sources. By heating molecular gas over large velocity widths, shocks can dramatically increase the efficiency of CO radiation. One source of shocks are outflows from young stars or star associations, examples in our Galaxy include the W 51, W 49 and DR 21 regions (Jaffe, Harris, & Genzel, 1987; Jaffe et al., 1989), in which bright CO $`\mathrm{J}=76`$ emission is observed with $``$20–50 $`\mathrm{km}\mathrm{s}^1`$ linewidths. In these sources, the CO luminosity in the broad-line shock-heated outflow components dominates that of the narrow-line PDR component. Another clue is the bright SiO observed in the central 150 pc of NGC 253 (Garcia-Burillo et al., 2000); the large implied abundance of SiO is attributed to outflows from massive young stars in a manner similar to the Galactic regions W 51 and DR 21. However, the comparison cannot apply globally in the nucleus of NGC 253. First, the excitation conditions in these Galactic shock regions are much higher that what we observe for the average properties of the nucleus of NGC 253. Observations of CO $`\mathrm{J}=21`$, $`\mathrm{J}=76`$ and $`\mathrm{J}=1615`$ toward W 51 IRS 2, for example, indicate that the gas temperature are 500–1000 K, unlikely for the bulk of the gas in NGC 253 given the rotational H2 observations suggesting lower temperatures. More importantly, in NGC 253, the near-IR H2 emission is characteristic of UV fluorescence, rather than the thermal spectrum produced by shock heating (Engelbracht et al., 1998) found in the Galactic outflow sources.

Shocks can also be generated with large-scale gas motions: supernova blast waves or spiral density waves. The nuclear starburst in NGC 253 is approximately bounded by an inner molecular ring with diameter 120 pc which corresponds to the inner Inner Linblad Resonance (iILR) of the galaxy’s bar. The large non-circular motions observed (as much as 40–50 $`\mathrm{km}\mathrm{s}^1`$) are attributed to gas accumulation and orbit crowding (see Buta & Combes (1996)), in the deep potential of the stellar mass, and it is plausible imagine that shocks from cloud collisions would be common. However, the difficulty again is that if the majority of the CO luminosity is to come from cloud collision shocks, they must be slow ($`v10`$ km s-1), non-dissociative shocks in order to reproduce the observed CO line ratios, and not excite the ro-vibrational H2 lines with a thermal spectrum. At these velocities, the collision timescale and energy available for collisions between clumps becomes small. Considering the molecular gas to come from assumed r$`2\mathrm{pc}`$, $`\mathrm{n}_{\mathrm{H}_2}=4.5\times 10^4`$ clumps, the mean separation of clumps $`\lambda `$ is around 25 pc. With a relative velocity of 10 $`\mathrm{km}\mathrm{s}^1`$ collisions will occur on timescales of

$$\tau \lambda ^3/(\sigma \mathrm{v})1\times 10^8\mathrm{yr}$$
(2)
with an available energy $`E=\frac{1}{2}M_{\mathrm{clump}}v^2`$. Such collisions, then, could in principle result in a luminosity to mass ratio of

$$\left(\frac{L}{M}\right)_{\mathrm{dyn}}=\frac{v^2}{2\tau }=\frac{v^3\sigma }{2\lambda ^3}8\times 10^5\frac{L_{}}{M_{}},$$
(3)
which is much smaller than the mass to luminosity relation implied by the model fit to the CO intensities ($`0.1`$). Lower densities and smaller clumps sizes mean more numerous clumps and a smaller mean separation. This can increase the luminosity by a factor of only 10–20, but it is still several orders of magnitude too small. Thus it is unlikely that a substantial fraction of the CO luminosity in the central $``$100 pc may be due to shocks from cloud-cloud collisions. Again, while higher-velocity shocks could provide much more energy, and the $``$40 $`\mathrm{km}\mathrm{s}^1`$ non-circular velocities are suggestive that such shocks could exist, the H2 ro-vibrational lines are inconsistent with any significant heating in fast shocks. Thus while some amount of shock heating is required to produce the high SiO abundance, it is not likely to be a dominant factor in heating the bulk of the molecular gas.

#### 4.1.3 Cosmic Ray Heating

It is widely believed that low energy cosmic-ray protons dominate the heating of the quiescent, UV-shielded molecular cloud cores in the disk of the Galaxy, keeping the gas temperature above that of the CMB (Goldsmith, & Langer, 1978, and references therein). A much larger cosmic ray density generated by massive star formation has also been invoked to explain the 50–150 K measured molecular gas temperatures in the central regions of M82 (SAH). Since the H2 cosmic ray dissociation cross sections are so small ($`3\times 10^{26}\mathrm{cm}^2`$), cosmic rays easily penetrate molecular cloud cores. The heating rate for cosmic rays can be characterized by a ionization rate $`\zeta _p`$ per H2 molecule, and an average $`\mathrm{\Delta }Q`$ deposited as thermal energy per ionization, that is:

$$\chi _{\mathrm{H}_2}=\zeta _p\mathrm{\Delta }Q.$$
(4)
Goldsmith, & Langer (1978) adopt values for $`\mathrm{\Delta }Q`$ between 17 eV and 20 eV. For the ionization rate we solve Eq. 3 of SAH which scales the cosmic ray density from that of the Galactic disk according to the supernova rate per unit area $`\psi `$ and a typical wind velocity $`v`$ by which cosmic rays are removed from the disk of a galaxy:

$$\frac{\zeta _{p,\mathrm{N253}}}{\zeta _{p\mathrm{Gal}}}=\frac{n_{cr\mathrm{N253}}}{n_{cr,\mathrm{Gal}}}\frac{\psi _{\mathrm{N253}}}{\psi _{\mathrm{Gal}}}\frac{v_{\mathrm{Galdisk}}}{v_{w,\mathrm{N253}}}.$$
(5)
The Galactic ionization rate $`\zeta _{p,\mathrm{Gal}}`$ is between $`2\times 10^{17}\mathrm{s}^1`$ (Goldsmith, & Langer, 1978) and $`7\times 10^{17}\mathrm{s}^1`$ (van Dishoeck, & Black, 1986), (see also Doty et al. (2002)). For the supernova rate in NGC 253 $`\psi _{\mathrm{N253}}`$, we take the lower limit derived for the central 360$`\times `$90 pc of 0.1 yr-1 (Antonucci & Ulvestad, 1988), giving $`3\mathrm{y}\mathrm{r}^1\mathrm{kpc}^2`$. To compare with the Galaxy, we take a rate of 5.8 per century over the entire disk (van den Bergh & Tammann, 1991, and references therein), or $`8\times 10^5\mathrm{yr}^1\mathrm{kpc}^2`$. The diffusion velocity in the Galactic disk is 10$`\mathrm{km}\mathrm{s}^1`$ (SAH), and the superwind velocity in NGC 253 is $``$500 $`\mathrm{km}\mathrm{s}^1`$, (Strickland et al., 2002). The ionization rate can then be scaled from that of the Galaxy:

$$\zeta _{\mathrm{p},\mathrm{N253}}=750\zeta _{\mathrm{p},\mathrm{Gal}}=1.55.3\times 10^{14}\mathrm{s}^1.$$
(6)
With a $`\mathrm{\Delta }Q`$ of 20 eV per ionization, the power deposited per H2 molecule is then 5–18 $`\times 10^{25}\mathrm{erg}\mathrm{s}^1`$. This estimate represents a lower limit to the cosmic ray heating, yet it is very close to the measured luminosity per H2 molecule in all of the CO lines, $`7\times 10^{25}\mathrm{erg}\mathrm{s}^1`$ with our assumed CO abundance.

To check that the increased CR flux does not destroy H2 and other molecular constituents, we compare our derived ionization rate with Eq. 10 of SAH. This requires that the ionization rate be less than the H2 formation rate:

$$\zeta _{\mathrm{p}\mathrm{N253}}<1.5\times 10^{18}\mathrm{n}\mathrm{T}^{1/2},$$
(7)
where $`\mathrm{n}`$ is the total hydrogen density. For our parameters, this requires $`\zeta _{\mathrm{p}\mathrm{N253}}<7.4\times 10^{13}\mathrm{s}^1`$, a factor of 50 from our lower limit of $`\zeta _{\mathrm{p}\mathrm{N253}}`$.

It is conceivable the increased CR density would dissociate the CO in the cloud core. Farquhar, Millar, & Herbst (1994) have constructed a chemical model of dark clouds to investigate the effects of cosmic rays ionization on the ion-neutral chemistry. For steady state models, with cloud density $`2\times 10^4\mathrm{cm}^3`$ and visual extinction $`A_V30`$ they find CO remains the dominant reservoir for carbon for the entire range of cosmic ray ionization rates investigated, $`\zeta =1.3\times 10^{19}8\times 10^{14}\mathrm{s}^1`$. The neutral carbon abundance grows sharply as $`\zeta `$ grows to $`2\times 10^{14}`$ approaching 25% of the CO abundance at the largest cosmic ray fluxes investigated. These results are consistent with the large $`\mathrm{C}^0/\mathrm{CO}`$ abundance ratios observed in the nuclei of both NGC 253 (Harrison et al., 1995), and M82 (Schilke et al., 1993). Since cosmic ray heating provides the necessary molecular gas luminosity and naturally accounts for the large observed neutral carbon abundance, while other heating sources are insufficient by large factors, we favor favor the cosmic rays as the dominant heating source for the bulk of the molecular gas in the nucleus of NGC 253. Our analysis is fairly general, and we note that similar arguments could be made for other starburst nuclei where the rate of supernovae explosions approaches 0.1 per year in a $``$100 pc region.

### 4.2 Cloud Size, Relation to Atomic Gas

The total gas mass in a region compared with the derived density determines a volume filling factor of clouds. For our modeled dense component, the result is $`\varphi _V=2.6\times 10^3`$. The physical size of the clouds, however, is not constrained by our observations. The existence of HCN in dense cores (Paglione et al., 1995) indicates that many clouds must be large enough to shield their interiors from UV, approximately $`\mathrm{A}_\mathrm{V}>10`$ for the cloud radius. This puts a lower limit on the cloud size according to Eq (1). If we assume that the photodissociated gas, as traced through its \[CII\] and \[OI\] line emission forms a surface layer on the molecular cloud interiors, then we can estimate the cloud sizes. As discussed in Section 4.1.1, such a layer would have a thickness such that $`\mathrm{A}_\mathrm{V}3`$, and the atomic to molecular gas mass ratio requires that the atomic surface constitute $`3\times 10^2`$ of the total cloud mass. For spherical clouds, this implies a cloud radius corresponding to $`\mathrm{A}_\mathrm{V}300`$. Converting to physical sizes, these arguments then give $`\mathrm{r}_{\mathrm{cloud}}2.0\mathrm{pc}`$, at our modeled density and solar dust abundance. In a plane-parallel geometry, appropriate for clouds illuminated from a single direction, the geometric argument would suggest a cloud thickness corresponding to $`\mathrm{A}_\mathrm{V}100`$, or 0.7 pc. Adopting the 2 pc spherical clouds for argument, each has a mass of $`1.1\times 10^5`$ M, and there about 330 such clouds in the central 180 pc region. For reference, we note that increasing the CO abundance does not change the arguments above regarding the individual cloud geometry, but would correspond to a smaller molecular gas mass, lower volume filling factor, and fewer number of clouds. Table 2 summarizes the parameters for our adopted model.

## 5 Conclusions

We have detected bright CO $`\mathrm{J}=76`$ emission from the nucleus of NGC 253. We combine our observations with lower-J measurements of both $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ and an LVG model to estimate the CO column density and physical conditions of the molecular gas in the starburst nucleus. The column density of CO we infer based on the transitions above $`\mathrm{J}=10`$ is slightly larger than but similar to the estimates made by HHR99 for the central 23′′. However, we find that the excitation conditions of this gas must be very high to reproduce the observed $`\mathrm{J}=65`$ and $`\mathrm{J}=76`$ intensities. A single component with $`\mathrm{T}120\mathrm{K}`$, and $`\mathrm{n}4.5\times 10^4\mathrm{cm}^3`$ provides the best fit to all the observed transitions above $`\mathrm{J}=10`$ for both 12CO and 13CO, and is consistent with the ISO molecular H2 observations. We adopt this as the favored model. While a smaller additional mass of cool gas is required to account for the 13CO $`\mathrm{J}=10`$ intensity, and more cool gas may be present, the CO column density of warm molecular gas cannot be much less than our derived values.

The large mass of warm gas is very difficult to explain with a photodissociation region scenario. The molecular gas mass is a factor of 10–30 greater than the atomic gas mass, while the two should be comparable if heated with stellar UV as in the Galactic PDR templates. Dynamical heating of molecular gas from shocks is also unlikely based on energetic considerations, and given the lack of evidence for significant heating in fast shocks. On the other hand, the very large cosmic ray flux in the starburst nucleus of NGC 253 readily provides the necessary energy to heat the gas, and a means to deliver energy throughout the volume of large, UV-shielded clouds. It is likely the dominant heating term in the molecular ISM energy budget. While some of the warm molecular gas that we observe is likely heated by the far-UV radiation and shocks, our results strongly support the conjecture by (Suchkov, Allen, & Heckman, 1993) that the bulk of molecular gas in starburst nuclei is held at high excitation by the enhanced cosmic ray flux.

## 6 Acknowledgments

We are very grateful to the staff of JCMT for their support and patience during the first SPIFI observing run. We especially thank the director of the JCMT, Ian Robson for his support, and Wayne Holland for his tireless efforts in SPIFI’s setup and first observations. We are also indebted to Cornell engineer Chuck Henderson for solutions to difficult problems at the telescope on the instrument’s first observing run. We thank an anonymous referee for several helpful suggestions with this manuscript. This work was supported by NASA grants NAGW-4503 and NAGW-3925, and NSF grants OPP-8920223, and OPP-0085812.