text:
# Random hopping fermions on bipartite lattices: Density of states, inverse participation ratios, and their correlations in a strong disorder regime

## I Introduction

During the last decade, we have acquired much wisdom of Anderson localization And ; AALR . Multifractal scaling dimensions HJKPS ; LFSG , for example, have been calculated exactly for a random Dirac fermion model MCW ; Ber , and a lot of novel universality classes for disordered systems have been discovered Gad ; Zir ; AltZir ; SFBN ; BCSZ ; ASZ . Among the universality classes found so far, chiral (or sublattice) classes Gad ; Zir  have been attracting much interest Gad ; Fur ; HWKM ; Fuk ; AltSim ; GLL ; FabCas , since they are exactly on random critical points. Gade has predicted for these classes that the density of states (DOS) diverges at the critical points Gad .

Recently, Motrunich et.al. MDH  have claimed that such DOS as was predicted within a naive weak disorder approach is incorrect if one considers strong disorder features of the model. They have predicted an alternative expression of the DOS, which shows a bit weaker divergence.

On the other hand, Carpentier and Le Doussal CarDouN ; CarDouE  have proposed a renormalization group (RG) method to analyze strong disorder features at low temperatures of the XY gauge glass model CarOst ; RSN ; XYCal ; OzeNis ; XYDev ; Sch . Introducing higher charge fugacities, they have derived generalized RG equations for the Coulomb gas model, which is remarkably intimate relationship with a nonlinear diffusion equation called KPP equation KPP ; Bra ; EbeSaa . Their idea has been applied to a Dirac fermion model with a random vector field, and properties of a strong disorder regime as well as a weak disorder regime are successfully described in a unified way HorDou ; FukD .

Interestingly, Mudry et. al. MRF  have shown that the same RG method can apply to the random hopping problem: It has turned out that the disorder strength flows to a strong disorder regime, which provides the same DOS as Motrunich et. al. have predicted.

Recent numerical calculations RyuHat  for a random hopping fermion model on a square lattice with $`\pi `$ flux suggest the divergent DOS, and we can now believe that the model actually belongs to the chiral orthogonal class. Therefore, we would like to confirm as the next step whether strong disorder features of the chiral orthogonal class reveal themselves. Unfortunately, the DOS may not be a good quantity to grasp them, since the divergent behavior itself is technically hard to observe in numerical computations.

In order to clarify the strong disorder features of the model furthermore, we calculate, in this paper, some alternative quantities, the inverse participation ratios (IPR) and their spatial correlations, where such features also appear manifestly.

This paper is organized as follows: In the next section II, we introduce lattice models and their field theory in the scaling limit. In Sec. III, we derive the scaling equations taking into account higher powers of energy terms. In Sec. IV, we convert the scaling equations into KPP-type equation and discuss its asymptotic solution. By the use of them, we then calculate the DOS, IPR, and their spatial correlations in Sec. V. We give concluding remarks in Sec. VI.

## II Model

In this section, we first introduce the following two kinds of lattice models with sublattice symmetry as well as time-reversal symmetry: One is the model studied by Hatsugai et. al. HWKM  defined on a square lattice with $`\pi `$ flux, and the other is the one on a honeycomb lattice Sem . Next we derive a Dirac fermion model in the continuum limit of the lattice models at the band center.

### II.1 Lattice models

We consider the following non-interacting random hopping fermions on bipartite lattices in two dimensions:

$`H={\displaystyle \underset{i,j}{}}t_{ij}c_i^{}c_j,`$ (1)
where $`i,j`$ stands for summation over nearest neighbor pairs and $`t_{ij}=t_{0,ij}+\delta t_{ij}`$ is pure and random hopping amplitude, both of which can be chosen as real numbers due to time-reversal symmetry. The model defined on a square lattice with $`\pi `$ flux is here described with a gauge $`t_{0,ij}=()^{j_y}t_0`$. In the case of the honeycomb lattice, $`t_{0,ij}`$ is uniform, just defined as $`t_{0,ij}=t_0`$.

We are interested in the ensemble-average of the correlation functions

$`P^{(q_1,q_2)}(ij)`$ $`=\overline{\mathrm{\Psi }_n|c_i^{}c_i|\mathrm{\Psi }_n^{q_1}\mathrm{\Psi }_n|c_j^{}c_j|\mathrm{\Psi }_n^{q_2}}`$ (2)
$`=\overline{|\mathrm{\Psi }_n(i)|^{2q_1}|\mathrm{\Psi }_n(j)|^{2q_2},}`$
where $`|\mathrm{\Psi }_n=_i\mathrm{\Psi }_n(i)|i`$ with $`|i=c_i^{}|0`$ is a normalized wave functions with an energy eigenvalue $`E_n`$. Especially $`P^{(q,0)}P^{(q)}`$ defines the moments of the wavefunction, the IPR. To calculate the above correlation functions, the following Green functions may be useful:

$`\mathrm{\Gamma }(i)=\text{Im}i|{\displaystyle \frac{1}{ti(\omega iE)}}|i,`$
$`\mathrm{\Gamma }^{(q_1,q_2)}(i,j)=\mathrm{\Gamma }^{q_1}(i)\mathrm{\Gamma }^{q_2}(j).`$ (3)
Actually we have

$`\omega ^{q_1+q_21}\mathrm{\Gamma }^{(q_1,q_2)}(i,j)\underset{(\omega +0)}{}`$
$`C_{q_1+q_2}{\displaystyle \underset{n}{}}|\mathrm{\Psi }_n(i)|^{2q_1}|\mathrm{\Psi }_n(j)|^{2q_2}\delta (EE_n),`$ (4)
where numerical constant $`C_q`$ is given by $`C_q=\pi (2q3)!!/(2n2)!!`$. Then, the averaged DOS is given by $`\rho (E)=\overline{\mathrm{\Gamma }^{(1,0)}}`$. In the following, we will take the scaling limit of the lattice model near the zero energy and calculate the above correlations by the use of the field theory. In that case, it is difficult to define normalized wave functions which is crucial to the IPR. To avoid this, we alternatively define

$`\stackrel{~}{P}^{(q_1,q_2)}(ij)={\displaystyle \frac{\omega ^{q_1+q_21}}{C_{q_1+q_2}}}{\displaystyle \frac{\overline{\mathrm{\Gamma }^{(q_1,q_2)}(i,j)}}{\rho }},`$ (5)
with the normalization condition relaxed. Using these, we can define generalized IPR as $`\stackrel{~}{P}^{(q)}\stackrel{~}{P}^{(q,0)}`$. Note that $`\rho `$ in the denominator in Eq. (5) is needed for the “normalization” of the generalized IPR, i.e., $`\stackrel{~}{P}^{(1)}=1`$. In what follows, we calculate Eq. (5) using continuum theory and the RG method.

### II.2 Continuum limit

In the continuum limit around the zero energy of the above lattice models, via $`aix`$ and $`c_i/\sqrt{a}_je^{ik_{\mathrm{F}j}}\psi _j(x)`$ with a lattice constant $`a`$, we have the following action HWKM ; Fuk ; GLL ; MRF , after the choice of a suitable basis

$`S={\displaystyle d^2x\overline{\psi }_i\left[i\gamma _\mu \left(_\mu +A_\mu \right)+iM_1+M_2\gamma _5\right]\psi _i},`$ (6)
where $`\gamma _1=\sigma _1`$, $`\gamma _2=\sigma _2`$, and $`\gamma _5=\sigma _3`$ are usual Pauli matrices, and $`i=1,2`$ denotes a “flavor” due to the species doubling. The probability distributions of the randomness are

$`P[A_\mu ]\mathrm{exp}\left({\displaystyle \frac{1}{2g_A}}{\displaystyle d^2xA_\mu ^2}\right),`$
$`P[M_\mu ]\mathrm{exp}\left({\displaystyle \frac{1}{2g_M}}{\displaystyle d^2xM_\mu ^2}\right).`$ (7)
The difference between two lattice models introduced above is that the initial strength $`g_A`$ and $`g_M`$ is given by $`g_A=g_M`$ for the model on a square lattice while $`g_Ag_M`$ for the one on a honeycomb lattice. However, it has nothing to do with the long distance behavior, since $`g_A`$ only is renormalized, as we shall see soon.

The energy term in Eq. (3) are

$`S_y=y_1{\displaystyle d^2x𝒴},`$ (8)
where $`y_1=\omega iE`$ and $`𝒴=𝒪_++𝒪_{}`$ with

$`𝒪_+=\overline{\psi }_{R1}\overline{\psi }_{L2}\overline{\psi }_{R2}\overline{\psi }_{L1}`$
$`𝒪_{}=\psi _{L2}\psi _{R1}\psi _{L1}\psi _{R2},`$ (9)
Here, chiral fermions have been defined by

$`\overline{\psi }_i=i\left(\begin{array}{cc}\overline{\psi }_{Li},\hfill & \overline{\psi }_{Ri}\hfill \end{array}\right),\psi _i=\left(\begin{array}{c}\psi _{Ri}\hfill \\ \psi _{Li}\hfill \end{array}\right).`$ (13)
In the following section, we add powers of this energy term to the action by the use of the replica method. They play a crucial role in the RG analysis, as we shall see.

## III Renormalization group analysis

In the previous work within usual weak disorder approaches, the replica method Gad ; Fuk ; FabCas  and supersymmetry method GLL  have been used to take quenched average. The latter method is quite interesting, since it enables us to derive exact beta functions for the present model. The recent work by Mudry et. al. has also used the supersymmetry method, but with higher energy (or “fugacity” in the context of the Coulomb gas model by Carpentier and Le Doussal) terms being added, anomalous dimensions of them may be hard to obtain beyond one-loop order.

Therefore, we will use the replica method for simplicity and derive the RG equations up to one-loop order for consistency. We follow a similar formulation developed by Mudry et. al. MRF .

### III.1 Scaling equations for fugacities

Let us introduce the replica, $`\psi _i\psi _{ia}`$ with a replica number $`m`$; $`a=1,2,\mathrm{},m`$. Then ensemble-averaged correlation functions can be calculated as

$`\overline{\mathrm{\Gamma }^{(q_1,q_2)}(xy)}`$
$`=𝒴_{a_1}𝒴_{a_2}\mathrm{}𝒴_{a_{q_1}}(x)𝒴_{b_1}𝒴_{b_2}\mathrm{}𝒴_{b_{q_2}}(y),`$ (14)
in the replica limit $`m0`$.

We utilize the Hodge decomposition to take ensemble-average over the vector field,

$`A_\mu =ϵ_{\mu \nu }_\nu \phi +_\mu \varphi ,`$ (15)
which converts the probability distribution (7) into

$`P[A]\mathrm{exp}\left\{{\displaystyle \frac{1}{2g_A}}{\displaystyle d^2x\left[\left(_\mu \phi \right)^2+\left(_\mu \varphi \right)^2\right]}\right\}.`$ (16)
The gauge transformation

$`\overline{\psi }_{ia}\overline{\psi }_{ia}e^{i\phi \gamma _5+\varphi },`$
$`\psi _{ia}e^{i\phi \gamma _5\varphi }\psi _{ia},`$ (17)
together with $`M_1\pm iM_2e^{\pm 2i\phi }(M_1\pm iM_2)`$, which gives rise to no change to the probability distribution $`P[M_\mu ]`$, yields the replicated action

$`S^{(m)}={\displaystyle d^2x\left(_0+_1+2g_M𝒪_M+\underset{n=1}{\overset{m}{}}y_n𝒴^n\right)},`$ (18)
with

$`_0={\displaystyle \frac{1}{2g_A}}\left(_\mu \varphi \right)^2,`$
$`_1=\overline{\psi }_{R\alpha }2_{\overline{z}}\psi _{R\alpha }+\overline{\psi }_{L\alpha }2_z\psi _{L\alpha },`$
$`𝒪_M=J_{R\alpha \beta }J_{L\beta \alpha },`$ (19)
where we have denoted $`\alpha =ia`$ and $`\beta =jb`$ for simplicity, and defined currents as

$`J_{R\alpha \beta }=\overline{\psi }_{R\alpha }\psi _{R\beta },J_{L\alpha \beta }=\overline{\psi }_{L\alpha }\psi _{L\beta }.`$ (20)
Note that the field $`\phi `$ has been decoupled and hence omitted in Eq. (18).

The last term in Eq. (18) is the $`n`$th power of the energy term, $`𝒴^n=_{\{a\}}𝒴_{a_1}\mathrm{}𝒴_{a_n}`$. Expanding these with respect to $`O_\pm `$ defined in Eq. (9), it turns out that the most relevant operators are $`_{\{a\}}\left(𝒪_{+a_1}\mathrm{}𝒪_{+a_n}+𝒪_{a_1}\mathrm{}𝒪_{a_n}\right)𝒪_+^n+𝒪_{}^n𝒴_{\mathrm{mr}}^n`$, which is actually converted via the gauge transformation (17) into

$`𝒴_{\mathrm{mr}}^n=e^{2n\varphi }𝒪_+^n+e^{2n\varphi }𝒪_{}^n.`$ (21)
As we shall see momentarily, it is important to include all these terms even if we started initially with $`𝒴_1`$ only, since these terms are generated in the process of the renormalization. To see this, note that the dimension of the vertex operator $`e^{n\varphi }`$ is $`n^2g_A/\pi `$ and therefore, the following OPE holds;

$`e^{n\varphi (z)}e^{n^{}\varphi (0)}{\displaystyle \frac{1}{|z|^{2nn^{}g_A/\pi }}}e^{(n+n^{})\varphi (0)}.`$ (22)
This equation is actually involved with a “fusion” of the energy terms, generating the higher powers of operators $`𝒴^n`$.

Based on the above OPE, let us now compute the OPE between $`𝒴_{\mathrm{mr}}^n`$ and $`𝒴_{\mathrm{mr}}^n^{}`$. One problem is the cross term such as $`e^{n\varphi (z)}O_+^n(z)e^{n^{}\varphi (0)}O_{}^n^{}(0)`$. Its most singular term in the OPE yields the exponent $`2nn^{}g_A/\pi +|nn^{}|`$, and therefore at least when $`g_A\pi /2`$, such OPE gives rise to no singularity. Therefore, we will neglect these terms, and then we reach the following OPE for $`𝒴_{\mathrm{mr}}^n`$,

$`𝒴_{\mathrm{mr}}^n(z)𝒴_{\mathrm{mr}}^n^{}(0){\displaystyle \frac{1}{|z|^{2nn^{}g_A/\pi }}}𝒴_{\mathrm{mr}}^{n+n^{}}(0).`$ (23)
This leads to the scaling equations

$`{\displaystyle \frac{dy_n}{dl}}=\beta _n(2x_n)y_n\pi {\displaystyle \underset{n^{}=1}{\overset{m}{}}}y_n^{}y_{nn^{}},`$ (24)
where $`x_n`$ is the scaling dimension of the operator $`𝒴_{\mathrm{mr}}^n`$ given by $`x_n=n\frac{g_A}{\pi }n^2`$ for the present. However, the $`𝒪_M`$ yields an anomalous dimension, which is determined in the next subsection. The correlation function $`\mathrm{\Gamma }^{(q_1,q_2)}(xy)`$ is dominated by

$`\overline{\mathrm{\Gamma }^{(q_1,q_2)}(xy)}=𝒴_{\mathrm{mr}}^{q_1}(x)𝒴_{\mathrm{mr}}^{q_2}(y).`$ (25)
In Sec. V, we will determine its long distance behavior using the RG equation for these correlation functions.

### III.2 Scaling equations for random mass and vector-field couplings

In this subsection, we derive the scaling equation for $`g_M`$. The basic OPE for the currents are GLL

$`J_{R\alpha \beta }(z)J_{R\kappa \lambda }(w){\displaystyle \frac{\delta _{\alpha \lambda }\delta _{\beta \kappa }}{4\pi ^2(zw)^2}}+{\displaystyle \frac{1}{2\pi (zw)}}\left[\delta _{\beta \kappa }J_{R\alpha \lambda }(w)\delta _{\alpha \lambda }J_{R\kappa \beta }(w)\right],`$
$`J_{R\alpha \beta }(z)\overline{\psi }_{R\kappa }{\displaystyle \frac{\delta _{\beta \kappa }}{2\pi (zw)}}\overline{\psi }_{R\alpha },J_{R\alpha \beta }(z)\psi _{R\kappa }{\displaystyle \frac{\delta _{\alpha \kappa }}{2\pi (zw)}}\psi _{R\beta },`$ (26)
and similar for $`L`$-movers. Therefore, we have the following OPE for $`𝒪_M`$ GLL

$`𝒪_M(z)𝒪_M(0){\displaystyle \frac{1}{4\pi ^2|z|^2}}\left[2𝒪_A(0)4m𝒪_M(0)\right],`$
(27)
$`𝒪_M(z)𝒪_A(0)𝒪_A(z)𝒪_A(0)0,`$ (28)
$`𝒪_M(z)𝒴_n(0){\displaystyle \frac{n}{4\pi ^2|z|^2}}\left[𝒴_n(0)+\mathrm{\Delta }𝒴_n(0)\right].`$ (29)
Here, another operator $`𝒪_AJ_{R\alpha \alpha }J_{L\beta \beta }`$ is generated. Note, however, that it is just an operator appearing after ensemble-average directly over $`A_\mu `$ instead of the gauge transformation (17), and hence it is involved in the renormalization of $`g_A`$. Note also that the OPE (29) contains operators not included in $`𝒴_n`$. However, we neglect them in order to close the OPE algebra. The RG equations for $`g_A`$ and $`g_M`$ in the replica limit $`(m0)`$ are then given by

$`{\displaystyle \frac{dg_A}{dl}}=\beta _A{\displaystyle \frac{g_M^2}{\pi }},`$
$`{\displaystyle \frac{dg_M}{dl}}=\beta _M0,`$ (30)
and the anomalous dimension of $`𝒴_{\mathrm{mr}}^n`$ is calculated as

$`x_n=\left(1{\displaystyle \frac{g_M}{\pi }}\right)n{\displaystyle \frac{g_A}{\pi }}n^2.`$ (31)
If we neglected the fusion of the energy terms, this would serve as the scaling dimension of the operator $`𝒴_{\mathrm{mr}}^n`$ and govern the long distance behavior of the correlation function $`𝒴^n`$.

## IV KPP equation

So far we have derived the scaling equations for $`y_n`$, $`g_A`$, and $`g_M`$. Those for $`g_A`$ and $`g_M`$ are actually scaling equations in the replica limit. Basically, we also have to solve the equations for $`y_n`$ and take the replica limit. To this end, we utilize the method developed by Carpentier and Le Doussal CarDouN .

Define the following distribution function $`P(l,u)`$ of $`y_n`$

$`y_n(l)`$ $`=`$ $`{\displaystyle \frac{2}{\pi n!}}{\displaystyle 𝑑ue^{nu}P(u,l)}`$ (32)
$``$ $`{\displaystyle \frac{2}{\pi n!}}e^{nu}_P.`$
Then it turns out that the resultant equation for $`P(u,l)`$ is free from $`m`$, and we need not to worry about the replica limit. Moreover, define

$`G(x,l)`$ $`=1\mathrm{exp}\left[e^{ux+(1g_M/\pi )l}\right]_P.`$ (33)
It is readily seen that this function $`G`$ obeys the following KPP-type equation,

$`{\displaystyle \frac{1}{2}}_lG=D(l)_x^2G+G(1G),`$ (34)
where $`D(l)=\frac{g_A(l)}{2\pi }`$ is a diffusion constant. The boundary condition is $`G(\mathrm{},0)=1`$ and $`G(\mathrm{},0)=0`$. The difference from the normal KPP equation is that the diffusion constant $`D`$ depends on $`l`$. In analyzing the above equation, we take an “adiabatic” approximation as Mudry et. al. applied to the same problem. Namely, we first neglect the $`l`$-dependence of the diffusion constant $`D`$ for a while, and after solving the KPP equation, we take the $`l`$-dependence into account. This scheme may hold at least in the leading approximation. It should be noted that the function $`G`$ is a generating functional of $`y_n`$. Actually, expanding with respect to $`e^x`$ leads to

$`G(x,l)={\displaystyle \frac{\pi }{2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}()^ny_n(l)e^{n[x(1g_M/\pi )l]}.`$ (35)
In the limit $`l\mathrm{}`$, the KPP equation asymptotically allows the following traveling wave solution $`G(x,l)G\left(xm(l)\right)`$, where $`m(l)`$ is a velocity (times $`l`$) of the traveling wave. In order to determine the dimension of the operator $`𝒴_{\mathrm{mr}}^q`$, let us perturb the action with it, setting the initial coupling $`y_n=y_q\delta _{n,q}`$ in Eq. (18). Then it turns out that

$`m(l)=\{\begin{array}{cc}2(Dq+1/q)l+O(1),\hfill & \text{for}q<1/\sqrt{D},\hfill \\ 4\sqrt{D}\left(l\frac{1}{8}\mathrm{ln}l\right)+O(1),\hfill & \text{for}q=1/\sqrt{D},\hfill \\ 4\sqrt{D}\left(l\frac{3}{8}\mathrm{ln}l\right)+O(1),\hfill & \text{for}q>1/\sqrt{D},\hfill \end{array}`$ (39)
The average value of $`u`$ can be evaluated as CarDouN ; FukD

$`u_Pm(l)\left(1{\displaystyle \frac{g_M}{\pi }}\right)l,`$ (40)
Using this, we define the typical value of $`y_n`$ as

$`y_{\mathrm{typ},n}(l)`$ $`={\displaystyle \frac{2}{\pi n!}}e^{nu_P}.`$ (41)
It is readily seen that $`y_{\mathrm{typ},n}`$ satisfies the following “scaling equation”,

$`{\displaystyle \frac{dy_{\mathrm{typ},n}}{dl}}`$ $`=n\left[{\displaystyle \frac{dm(l)}{dl}}\left(1{\displaystyle \frac{g_M}{\pi }}\right)\right]y_{\mathrm{typ},n}`$ (42)
$`=z_ny_{\mathrm{typ},n},`$
To be explicit, Eq. (39) yields

$`z_n`$ $`=\{\begin{array}{cc}2\left(1\frac{g_M}{\pi }\right)n+2Dn^2,\hfill & \text{for}n<1/\sqrt{D},\hfill \\ 4\sqrt{D}\left(1\frac{1}{8l}\right)n\left(1\frac{g_M}{\pi }\right)n,\hfill & \text{for}n=1/\sqrt{D},\hfill \\ 4\sqrt{D}\left(1\frac{3}{8l}\right)n\left(1\frac{g_M}{\pi }\right)n,\hfill & \text{for}n>1/\sqrt{D}.\hfill \end{array}`$ (46)
Notice that the first line of the equation reproduces the naive dynamical exponent of $`y_n`$, i.e., $`z_n=2x_n`$, where $`x_n`$ is given by Eq. (31). We denote the exponent thus defined as $`z_{\mathrm{w},n}`$, which we shall discuss in Sec. V.3. This result implies that effective scaling dimensions of $`y_n`$ can be modified due to the fusion of the energy terms.

Let us now take account of the dependence of $`D`$ on $`l`$ through $`g_A`$. It should be stressed that even if one takes into account the $`l`$-dependence of $`D`$ in deriving Eq. (46) from Eq. (42), the following calculations hold. For sufficiently large $`l`$, $`z_n`$ is given by the last of the selection rule above, since $`1/\sqrt{D}1/\sqrt{l}1`$. Namely, we have

$`z_n`$ $``$ $`4n\sqrt{D}`$
$``$ $`nz(l),`$
$`z(l)`$ $``$ $`\text{const.}\left({\displaystyle \frac{g_M}{\pi }}\right)l^{\frac{1}{2}},`$ (47)
and using this, we define the typical scaling exponent $`x_{\mathrm{typ},n}`$ of the operator $`𝒴_{\mathrm{mr}}^n`$ as

$`x_{\mathrm{typ},n}(l)=2z_n2nz(l),(l\mathrm{}).`$ (48)
There are two differences between $`z_n`$ and $`z_{\mathrm{w},n}`$: One is the $`l`$-dependence, which is essential to the behavior of the DOS, and the other is the $`n`$-dependence, which plays a role in the spatial correlation functions of the IPR.

## V The RG equations

In order to calculate the DOS, IPR, and their spatial correlations, we utilize the RG equations for the correlation functions $`\overline{\mathrm{\Gamma }^{(q_1,q_2)}(xy)}`$. Unfortunately, the equations are too hard to solve, so that we calculate their “typical” values based on the scaling equations for typical exponents derived in the previous section.

### V.1 One-point functions

In this subsection, we first derive the RG equation for $`\overline{\mathrm{\Gamma }^{(q,0)}}`$ to calculate the DOS and the IPR. Recall that $`𝒴_{\mathrm{mr}}^q`$ includes two fermion correlation functions such as $`\overline{\psi }_{R1a}(x)\overline{\psi }_{L2a}(x)`$, which diverges for a finite $`\omega `$. In order to get finite renormalized correlation functions, point-splitting the fields is necessary, $`\overline{\psi }_{R1a}(x)\overline{\psi }_{L2a}(x+a)`$. Then $`𝒴_{\mathrm{mr}}^q(x)`$ does not depend on $`x`$ but does on the cut-off $`a`$. Therefore under the scale transformation $`a(1+dl)a`$, it turns out that $`𝒴_{\mathrm{mr}}^q`$ obeys the RG equations

$`\left({\displaystyle \frac{}{l}}\beta _n{\displaystyle \frac{}{y_n}}\beta _A{\displaystyle \frac{}{g_A}}+\widehat{\gamma }_q\right)𝒴_{\mathrm{mr}}^q=0,`$ (49)
where the beta functions are defined in Eqs. (24) and (30), and the (matrix of) anomalous dimension $`\widehat{\gamma }`$ is given by

$`\widehat{\gamma }_q𝒴_{\mathrm{mr}}^q=x_q𝒴_{\mathrm{mr}}^q+2\pi {\displaystyle \underset{n=1}{\overset{mq}{}}}𝒴_{\mathrm{mr}}^{q+n},`$ (50)
with $`x_q`$ defined in Eq. (31) Ber ; Zam . Since these equations are coupled together, it may be difficult to solve and take the replica limit.

One of possible approximations is to neglect the fusion of the operators $`𝒴_{\mathrm{mr}}^q`$ in Eqs. (23) and (24). This approximation corresponds to the conventional weak disorder approaches. The anomalous dimension $`\widehat{\gamma }_q`$ becomes in this case trivially diagonal and the RG equation is readily integrated. The results thus obtained are summarized in Sec. V.3. However, this scheme is not able to grasp essential strong disorder features the present model should have.

An alternative way is to compute the “typical” values of the correlation function, since we have been able to calculate $`y_{\mathrm{typ},n}`$ in the last section. This includes strong disorder effects, as it should be. Assuming that these typical values $`y_{\mathrm{typ},n}`$ correspond to the typical values of the correlation functions, which we denote as $`𝒴_{\mathrm{mr}}^n_{\mathrm{typ}}`$, one can simplify the RG equations above for these typical values, since the scaling equations for $`y_{\mathrm{typ},n}`$ in Eq. (42) are “diagonalized”, and so is the matrix of anomalous dimension $`\widehat{\gamma }_n`$. The RG equations for the typical values should be

$`\left({\displaystyle \frac{}{l}}\beta _A{\displaystyle \frac{}{g_A}}+x_{\mathrm{typ},q}\right)𝒴_{\mathrm{mr}}^q_{\mathrm{typ}}=0.`$ (51)
We are then led to

$`𝒴_{\mathrm{mr}}^q_{\mathrm{typ}}\mathrm{exp}\left[{\displaystyle ^l}𝑑l^{}x_{\mathrm{typ},q}(l^{})\right].`$ (52)
On the other hand, recall that the energy $`E`$ is related to $`y_1`$ as $`y_1=\omega iE`$. Corresponding to $`y_{\mathrm{typ},1}`$, let us define the typical energy $`y_{\mathrm{typ},1}=\omega _{\mathrm{typ}}iE_{\mathrm{typ}}`$. From Eq. (42), it follows that

$`{\displaystyle \frac{\mathrm{\Lambda }_{\mathrm{typ}}}{E_{\mathrm{typ}}}}`$ $`\mathrm{exp}\left[{\displaystyle ^l}𝑑l^{}z(l^{})\right]`$ (53)
$`=\mathrm{exp}\left[\text{const.}\left({\displaystyle \frac{g_M}{\pi }}\right)l^{\frac{3}{2}}\right],`$
where $`\mathrm{\Lambda }_{\mathrm{typ}}`$ is a renormalized energy and const. denotes a nonuniversal positive constant. Hence, as a function of $`E_{\mathrm{typ}}`$, we can rewrite the typical DOS $`\rho _{\mathrm{typ}}=𝒴_{\mathrm{mr}}^1_{\mathrm{typ}}`$ into

$`\rho _{\mathrm{typ}}(E_{\mathrm{typ}}){\displaystyle \frac{\mathrm{\Lambda }_{\mathrm{typ}}}{E_{\mathrm{typ}}}}\mathrm{exp}\left[c\left(\mathrm{ln}{\displaystyle \frac{\mathrm{\Lambda }_{\mathrm{typ}}}{E_{\mathrm{typ}}}}\right)^\kappa \right],`$ (54)
with

$`c\text{const.}\left({\displaystyle \frac{\pi }{g_M}}\right),\kappa ={\displaystyle \frac{2}{3}}.`$ (55)
This has been derived for the first time by Motrunich et. al. MDH  with the help of a strong coupling expansion, and rederived by Mudry et. al. MRF  by using the supersymmetry method. What is responsible for this value of $`\kappa `$ is the velocity selection rule (39) or (46) and resultant dynamical exponents (47). Recent numerical calculations RyuHat  support, to be sure, the divergent behavior of the DOS, but the nonuniversal constant in Eq. (55) may be too small for us to observe the exponent $`\kappa `$ in the numerical calculations, unfortunately.

To overcome the difficulty, we calculate the IPR for general $`q`$, which could make the strong disorder features manifest. For general $`q`$, we have

$`\omega _{\mathrm{typ}}^{q1}𝒴_{\mathrm{mr}}^q_{\mathrm{typ}}\rho _{\mathrm{typ}}(l),`$ (56)
which follows from

$`(q1)z+x_{\mathrm{typ},q}`$ $`=`$ $`(q1)z+\left(2qz\right)`$ (57)
$`=`$ $`2z`$
$`=`$ $`x_{\mathrm{typ},1}.`$
Therefore, we reach the following IPR:

$`\stackrel{~}{P}_{\mathrm{typ}}^{(q)}(E)\text{const.}`$ (58)
This implies that the multifractal scaling exponent $`\tau (q)`$, defined by $`\stackrel{~}{P}^{(q)}L^{\tau (q)}`$, is given by

$`\tau (q)=0,`$ (59)
for the typical values. This result is also attributed to the velocity selection rule (46) and therefore the strong disorder features of the present model. For references, we summarize in Sec. V.3, the results of a naive weak disorder approach.

### V.2 Two-point functions

Constant IPR obtained in the last subsection usually implies a localized wavefunction. However, it is believed in general that the zero-energy wavefunction of chiral models is extended. To clarify this point, we calculate the spatial correlations of the IPR in this subsection. The correlation function $`\overline{\mathrm{\Gamma }^{(q_1,q_2)}(xy)}=𝒴_{\mathrm{mr}}^{q_1}(x)𝒴_{\mathrm{mr}}^{q_2}(y)`$ obeys the following OPE

$`𝒴_{\mathrm{mr}}^{q_1}(x)𝒴_{\mathrm{mr}}^{q_2}(y)C_{q_1q_2}(xy)𝒴_{\mathrm{mr}}^{q_1+q_2}(y).`$ (60)
This may be valid for $`a|xy|L`$, where $`L`$ is a system size, and we will determine the $`r|xy|`$-dependence of the coefficient $`C_{q_1q_2}`$. It is easily verified that the coefficient $`C_{q_1q_2}`$ obeys the RG equations,

$`\left[\left(r{\displaystyle \frac{}{r}}\beta _n{\displaystyle \frac{}{y_n}}\beta _A{\displaystyle \frac{}{g_A}}\right)C_{q_1q_2}(r)C_{q_1q_2}(r)\widehat{\gamma }_{q_1+q_2}\right]𝒴_{\mathrm{mr}}^{q_1+q_2}+(\widehat{\gamma }_{q_1}+\widehat{\gamma }_{q_2})𝒴_{\mathrm{mr}}^{q_1}𝒴_{\mathrm{mr}}^{q_2}=0,`$ (61)
where $`\widehat{\gamma }`$ is defined in Eq. (50). These coupled set of equations become tractable for the typical values:

$`\left(r{\displaystyle \frac{}{r}}\beta _A{\displaystyle \frac{}{g_A}}+X_{\mathrm{typ},q_1,q_2}\right)C_{\mathrm{typ},q_1q_2}(r)=0,`$ (62)
where $`X_{\mathrm{typ},q_1,q_2}=x_{\mathrm{typ},q1}+x_{\mathrm{typ},q_2}x_{\mathrm{typ},q_1+q_2}`$. For sufficiently large $`l`$, it follows from Eqs. (47) and (48) that $`X_{\mathrm{typ},q_1,q_2}2`$. Thus we obtain

$`\stackrel{~}{P}_{\mathrm{typ}}^{(q_1,q_2)}(xy){\displaystyle \frac{\text{const.}}{|xy|^2}},`$ (63)
where we have used that previous result (58). Remarkably, the result shows that the exponent is a constant, 2, independent of $`q_1`$ and $`q_2`$ as well as the disorder strength. Since the correlation functions show power-law behavior, the wavefunction should not be localized.

### V.3 Weak disorder approximation

So far we have obtained several formulas taking into account strong disorder effects. As has been discussed by Motrunich et. al. as well as Mudry et. al., the exponent $`\kappa `$ in the DOS is $`\kappa =2/3`$, which is different from $`\kappa =1/2`$ obtained previously by Gade with the help of a naive weak disorder approach Gad . In order to clarify the strong disorder features of the IPR and their spatial correlations obtained in the previous subsection, we calculate them within the conventional weak disorder approach, for reference.

Without considering the fusion process of higher energy terms, the dynamical exponents of $`𝒴^n`$ is given just by the first line of Eq. (46). Namely,

$`z_{\mathrm{w},n}2Dn^2n^2z_\mathrm{w}(l),`$
$`z_\mathrm{w}(l)\text{const.}\left({\displaystyle \frac{g_M}{\pi }}\right)l.`$ (64)
Here, the label w means exponents within a weak disorder approximation. Note the difference of the $`l`$-dependence between Eqs. (47) and (64): This reflects the exponent $`\kappa `$ in the DOS. Neglecting the fusion of the operators $`𝒴_{\mathrm{mr}}^n`$, we find that the RG equation for $`𝒴_{\mathrm{mr}}^q`$ is diagonal, given by Eq. (51) but with the naive scaling dimension $`x_q2q^2z_\mathrm{w}(l)`$ in Eq. (31). Therefore, together with

$`{\displaystyle \frac{\mathrm{\Lambda }}{E}}`$ $`\mathrm{exp}\left[{\displaystyle ^l}𝑑l^{}z_\mathrm{w}(l^{})\right]`$ (65)
$`=\mathrm{exp}\left[\text{const.}\left({\displaystyle \frac{g_M}{\pi }}\right)l^2\right],`$
we have the same DOS (54) but with

$`\kappa _\mathrm{w}={\displaystyle \frac{1}{2}}.`$ (66)
Next, let us consider the correlation functions for general $`q`$. The solution of the RG equations yields

$`\omega ^{q1}𝒴_{\mathrm{mr}}^q/\rho `$
$`\mathrm{exp}\left\{{\displaystyle ^l}𝑑l^{}\left[(q1)z_\mathrm{w}x_q+x_1\right]\right\}`$
$`=\mathrm{exp}\left[q(q1){\displaystyle ^l}𝑑l^{}z_\mathrm{w}\right],`$ (67)
from which, together with Eq. (65), we reach

$`\stackrel{~}{P}_\mathrm{w}^{(q)}(E)\left({\displaystyle \frac{E}{\mathrm{\Lambda }}}\right)^{q(q1)}.`$ (68)
Negative exponent in the above equation suggests that the present weak disorder approach is incorrect and the result in the last subsection (58) is more reasonable. Since we are dealing with the same model in the same formalism, this result also implies that the exponent $`\kappa `$ in the DOS should be $`\kappa =2/3`$ rather than $`1/2`$.

In a similar way, we reach,

$`\stackrel{~}{P}_\mathrm{w}^{(q_1,q_2)}(xy)`$ $`\mathrm{exp}\left[cq_1q_2\left(\mathrm{ln}|xy|\right)^2\right]`$ (69)
$`\times \stackrel{~}{P}_\mathrm{w}^{(q_1+q_2)}(E),`$
where $`c`$ is given by the same disorder-dependent constant as in Eq. (55).

## VI Concluding remarks

In this paper we have reconsidered random hopping fermion models, based on the recent developments of the RG method for a strong disorder regime. We have used the replica method and calculated the DOS, IPR, and their spatial correlations. We have been able to reproduce the DOS predicted by Motrunich et. al. as well as Mudry et. al. Moreover, we have calculated the IPR and their spatial correlations, whose behavior near the zero energy is also different from the ones predicted by using conventional weak disorder approach.

In the present calculations, including higher powers of the energy perturbation (or higher fugacity terms in the context of Coulomb gas model) plays a crucial role, which has been invented for the first time by Carpentier and Le Doussal and applied to the present model by Mudry et. al.

Although recent numerical calculations have supported the divergent DOS toward the zero energy, it may be difficult to distinguish the exponent $`\kappa `$. We hope that IPR and their correlations calculated in the present paper may be useful to establish the strong disorder features of the random hopping models.

###### Acknowledgements.

We would like to thank H. Suzuki and S. Ryu for valuable discussions.


format: text

text:
# Fine Structures of Shock of SN 1006 with the Chandra Observation

## 1 Introduction

Since the discovery of cosmic rays (Hess, 1912), the origin and acceleration mechanism up to $`10^{15.5}`$ eV (the “knee” energy) have been long-standing problems. A breakthrough came from the X-ray studies of SN 1006; Koyama et al. (1995) discovered synchrotron X-rays from the shells of this SNR, indicating the existence of extremely high energy electrons up to the knee energy produced by the first order Fermi acceleration. Further, Tanimori et al. (1998) confirmed the presence of high energy electrons with the detection of the TeV $`\gamma `$-rays, which are cosmic microwave photons up-scattered by high energy electrons (the inverse Compton process) in the north east shell (the NE shell) of SN 1006. The combined analysis of the synchrotron X-rays and inverse Compton TeV $`\gamma `$-rays nicely reproduces the observed flux and spectra, and predicts a rather weak magnetic field of 4–6 $`\mu `$G (Tanimori et al., 1998, 2001).

Since these discoveries, detection of synchrotron X-rays and/or TeV $`\gamma `$-rays, from other shell-like SNRs has been accumulating: G347.3$``$0.5 (Koyama et al., 1997; Slane et al., 1999; Muraishi et al., 2000; Enomoto et al., 2002), RCW 86 (Bamba et al., 2000; Borkowski et al., 2001b), and G266.6$``$1.2 (Slane et al., 2001). These discoveries provide good evidence for the cosmic ray acceleration at the shocked shell of SNRs. The mechanism of the cosmic ray acceleration has also been studied for a long time and the most plausible process is a diffusive shock acceleration (DSA) (Bell, 1978; Blandford & Ostriker, 1978; Drury, 1983; Blandford & Eichler, 1987; Jones & Ellison, 1991; Malkov & Drury, 2001).

Apart from the globally successful picture of DSA, detailed but important processes, such as the injection, magnetic field configuration, and the reflection of accelerated particles, have not yet been well understood. The spatial distribution of accelerated particles responsible for the non-thermal X-rays, may provide key information on these unclear subjects. Previous observations, however, are limited in spatial resolution for a detailed study on the structure of shock acceleration process and injection efficiency. Although many observations and theoretical models are made for SN 1006, these problems are still open issue (Reynolds, 1998; Aharonian & Atoyan, 1999; Vink et al., 2000; Ellison, Berezhko, & Baring, 2000; Dyer et al., 2001; Allen, Petre, & Gotthelf, 2001; Berezhko et al., 2002).

In this paper, we report on the first results of the spectral and spatial studies on the thermal and non-thermal shock structure in the NE shell of SN 1006 with Chandra (§ 3). In § 4.1 and § 4.2, we discuss the spectral analyses and determine the scale widths of the structures for thermal and non-thermal electrons on the base of a simple DSA with shock parallel magnetic field. We also derive the injection efficiency ($`\eta `$) of non-thermal electrons from the thermal plasma near the shock front. Based on these results, we discuss possible implications on the DSA process in the NE shell of SN 1006. In this paper, we assume the distance of SN 1006 to be 1.8 kpc (Green, 2001).

## 2 Observation

We used the Chandra archival data of the ACIS on the NE shell of SN 1006 (Observation ID = 00732) observed on July 10–11, 2000 with the targeted position at (RA, DEC) = ($`15^\mathrm{h}03^\mathrm{m}51\stackrel{\mathrm{s}}{\mathrm{.}}6`$, $`41^\mathrm{d}51^\mathrm{m}18\stackrel{\mathrm{s}}{\mathrm{.}}8`$). The satellite and instrument are described by Weisskopf et al. (1996) and Garmire et al. (2000), respectively. CCD chips I2, I3, S1, S2, S3, and S4 were used with the pointing center on S3. Data acquisition from the ACIS was made in the Timed-Exposure Faint mode with the readout time of 3.24 s. The data reductions and analyses were made using the Chandra Interactive Analysis of Observations (CIAO) software version 2.2.1. Using the Level 2 processed events provided by the pipeline processing at the Chandra X-ray Center, we selected ASCA grades 0, 2, 3, 4, and 6, as the X-ray events. High energy electrons due to charged particles and hot and flickering pixels were removed. The effective exposure was $`68`$ ks for the observation. In this paper, we concentrated on the data of S3 (BI chip) because this chip has the best efficiency in soft X-rays required for the spectral analyses and its on-axis position provides the best point-spread function required for the spatial analysis.

## 3 Analyses and Results

### 3.1 Overall Image

Figure 1 shows the true-color image from I2, I3, and S1–S4 for the NE shell of SN 1006. The image is contrasted in the 0.5–2.0 keV band (hereinafter, the soft1 band; red) and in the 2.0–10.0 keV band (hereinafter, the hard band; blue) and binned to a resolution of 1 arcsec. The fine spatial resolution of Chandra unveils extremely narrow filaments in the hard band. They are running from north to south along the outer edge of the NE shell, parallel to the shock fronts observed by H$`\alpha `$ emission line (Winkler & Long, 1997). These filaments resemble the sheet-like structure of the shock simulated by Hester (1987). The soft1 band image, on the other hand, has a larger scale width similar to the ROSAT HRI image (Winkler & Long, 1997). Many clumpy sub-structures are also seen in this energy band.

### 3.2 Inner Shell Region

To resolve the thermal and non-thermal components, we made a spectrum from a bright clump found in the soft1 band image, which is located in the inner part of the NE shell (“Inner region” with the dashed ellipse in Figure 2). The background region was selected from a region out of the SNR, as is shown in Figure 2 with the dashed lines.

The background-subtracted spectrum shown in Figure 3 has many emission lines. We hence determined the peak energies of the 5 brightest lines with a phenomenological model, a power-law continuum plus Gaussian lines. The most strong line structures are the peak at 0.55 keV and the hump at 0.67 keV. These energies are nearly equal to the K$`\alpha `$ and Lyman $`\alpha `$ lines of He- and H-like oxygen, hence are attributable to highly ionized oxygen. Likewise, the other clear peaks at 0.87, 1.31, and 1.76 keV are most likely He-like K$`\alpha `$ of Ne, Mg and Si, respectively. However, in detail all the observed line energies are systematically smaller than those of the relevant atomic data. These apparent energy shifts have been usually observed in a young SNR plasma in non-equilibrium ionization (NEI). The “energy shift” in this case is due to the different line ratio of many sub-levels and/or different ionization states. The oxygen Lyman $`\alpha `$ is isolated from the other lines of different ionization states, hence the NEI effect gives no energy shift. Still we see apparent down-shift of the observed line energy from that of the laboratory data. He-like K$`\alpha `$ lines are complex of many fine structures with the split-energy of at most $``$ 25 eV (for He-like silicon). Although the energy shift of He-like K$`\alpha `$ lines due to NEI should be smaller than this split-energy, the observed energy shifts are systematically larger than the split-energy. We therefore regard that the apparent energy shifts are due mainly to energy calibration errors, hence fine-tuned the energy gain to reduce by 3.8%, the average shift of the 5 brightest lines. We then fitted the spectrum with a thin thermal plasma model in NEI calculated by Borkowski et al. (2001a). The abundances of C, N, O, Ne, Mg, Si, S, and Fe in the plasma were treated to be free parameters, whereas those of the other elements were fixed to the solar values (Anders & Grevesse, 1989). The absorption column was calculated using the cross sections by Morrison & McCammon (1983) with the solar abundances. Since this NEI model exhibited systematic data excess at high energy above 2 keV, we added a power-law component and the fit improved dramatically. Figure 3 and Table 1 show the best-fit models (dashed and solid lines for thermal and power-law components) and parameters, respectively.

Instead of the phenomenological power-law model, we applied srcut in the XSPEC package as a more physical model. Details of the srcut model fitting are given in § A.1. The best-fit $`\nu _{\mathrm{rolloff}}`$ is 9.2 (8.6–10.3) $`\times 10^{16}`$ Hz, with better reduced $`\chi ^2`$ of 389.0/215 than that of the power-law model of 447.9/215 (see Table 1).

As for the thermal components, we also tried the fitting with a plane shock model (XSPEC model $`vpshock`$) plus either a power-law or srcut and found no essential difference from the case of an NEI model.

Although these simple models globally follow the data very well, all are rejected in the statistical point of view, leaving wavy residuals near the line structure as shown in Figure 3 (lower panel). This may be caused by improper response function in energy scale and/or in energy resolution. We assumed that the photons are uniformly distributed in flux and in temperature through the whole source region. This simple assumption may also partly responsible for the above systematic error, because in reality the source region is apparently clumpy (see Figure 1) and may have different temperatures, abundances, and/or ionization time scales. Since our principal aim of this paper is to examine the spatial structures and the spectra of non-thermal component, we do not examine in further detail for the thermal model. In the following analyses and discussion, we use the physical parameters cited in Table 1 as a good approximation.

The spectrum of the inner region clump is softer than any other regions in the NE shell, which indicates that the contribution of the thermal component is the largest. Nevertheless the thermal photons are only 0.02% of the non-thermal ones if we limit the energy band to 2.0–10.0 keV (the hard band); the non-thermal photons in the hard band are 8.9$`\times 10^2`$ cnts s-1, while thermal photons are 2.0$`\times 10^5`$ cnts s-1. Therefore, in the following spatial analyses, we regard that all the photons in the hard band are non-thermal origin.

As for the spatial analysis of the thermal emission, we use the limited band of 0.4–0.8 keV (hereafter; the soft2 band) to optimize the signal-to-noise ratio, in which K-shell lines from He-like oxygen (0.57 keV) contribute a most fraction of the X-ray emission (see Figure 3). Even in this thermal-optimized band, however, the count rates of the thermal and non-thermal emissions are comparable: thermal photons are 8.1$`\times 10^1`$ cnts s-1, while those of non-thermal are 5.4$`\times 10^1`$ cnts s-1.

### 3.3 The Filaments

The outer edge of the NE shell is outlined by several thin X-ray filaments. For the study of these filaments, we selected 6 rectangle regions in Figure 1, in which the filaments are straight and free from other structures like another filament and/or clumps. These regions (solid boxes) are shown in Figure 2 with the designations of No.1–6 from north to south. Since the SNR shell is moving (expanding) from the right to the left, we call the right and left side as downstream and upstream following the terminology of the shock phenomena.

Figure 4 shows the intensity profile in the hard (2.0–10.0 keV: upper panel) and soft2 (0.4–0.8 keV: lower panel) bands for each filament with the spatial resolution of 0.5 arcsec, where the horizontal axis ($`x`$-coordinate) runs from the east to west (upstream to downstream) along the line normal to the filaments. We see very fast decay in the downstream side and even faster rise in the upstream side.

To estimate the scale width, we define a simple empirical model as a function of position ($`x`$) for the profiles;

$$f(x)=\{\begin{array}{ccc}\hfill A\mathrm{exp}|\frac{x_0x}{w_\mathrm{u}}|& \mathrm{in}\mathrm{upstream}\hfill & \\ \hfill A\mathrm{exp}|\frac{x_0x}{w_\mathrm{d}}|& \mathrm{in}\mathrm{downstream},\hfill & \end{array}$$
(1)
where $`A`$ and $`x_0`$ are the flux and position at the filament peak, respectively. The scale widths are given by $`w_\mathrm{u}`$ and $`w_\mathrm{d}`$ for upstream and downstream, respectively (hereafter, “u” and “d” represent upstream and downstream, respectively). Since the scale width of the filaments is larger than the spatial resolution of Chandra ($`0.5`$ arcsec = 1 bin in Figure 4), we ignore the effect of the point-spread function.

#### 3.3.1 Non-Thermal Structure

As already noted in § 3.2, the hard band (2.0–10.0 keV) flux is nearly pure non-thermal origin (see also the next paragraph). We therefore used the hard band profiles for the study of the non-thermal X-ray structures. The hard band profiles were fitted with a function $`f^\mathrm{h}(x)+C^\mathrm{h}`$ (hereafter, “h” represents the hard X-ray band), where $`C^\mathrm{h}`$ is the background constant, which includes the cosmic and Galactic X-ray background and non X-ray events. The fittings were statistically accepted for all the filament profiles. The best-fit models and parameters are shown in Figure 4 with the solid lines and in Table 2, respectively.

We then made the spectra of the filaments within the scale widths: ($`x_0^\mathrm{h}w_\mathrm{u}^\mathrm{h}xx_0^\mathrm{h}+w_\mathrm{d}^\mathrm{h}`$) in Figure 4. The background spectra were made from the off-filament downstream regions. All the background-subtracted spectra are featureless (no line structure) and extend to high energy side, which were fitted with an absorbed power-law model with the best-fit parameters given in Table 3. We thus confirm that the hard X-ray profiles represent those of non-thermal X-rays.

To increase statistics, we summed all the data of 6 filaments (the combined-filament). The best-fit power-law model parameters for the spectrum of the combined-filament are listed in Table 3. The spectrum was also fitted with a srcut model (see § A.1). The best-fit $`\nu _{\mathrm{rolloff}}`$ and the other parameters are also listed in Table 3. We note that the srcut model give slightly better $`\chi ^2`$/d.o.f. than that of a phenomenological power-law model (see Table 3).

We further spatially divided the combined-filament and analyzed each spectrum. We however found no significant difference between the downstream and upstream, nor within the downstream; the photon index was nearly constant along the $`x`$-axis in the combined-filament.

#### 3.3.2 Thermal Structure

To examine the structure of the thermal components, we used the soft2 band profiles (the lower panels of Figure 4). Contamination from the non-thermal photons would be very large even in this optimized band (see § 3.2). Therefore we calculated the flux ratio of the non-thermal photons between in the hard and soft2 bands using the spectral parameters given in Table 3. From the flux ratio and the best-fit hard band profiles, we estimated the non-thermal contaminations as are given with the dashed lines in Figure 4.

After subtracting these non-thermal contaminations, we fitted the soft2 band profiles with a model of $`f^\mathrm{s}(x)+C^\mathrm{s}`$, where $`C^\mathrm{s}`$ is the background constant in the same sense as $`C^\mathrm{h}`$. Note that although we use “s” as the abbreviation of the soft2 band, it actually represents thermal X-rays.

From Figure 4, we see that most of the photons at the filament peak is non-thermal origin, hence the statistics becomes poor to determine the position of the thermal peak ($`x_0^\mathrm{s}`$) independently. We thus fixed $`x_0^\mathrm{s}`$ to the best fit peak in the hard hand $`x_0^\mathrm{h}`$. Also we set the upper bound of the fitting parameter $`w_\mathrm{d}^\mathrm{s}`$ to be 900 arcsec, the same as the radius of SN 1006 (Green, 2001). The best-fit models and parameters for the thermal components ($`f^\mathrm{s}(x)`$) are shown with dotted lines in the lower panels of Figure 4 and in Table 2, respectively.

#### 3.3.3 Non-Thermal versus Thermal

Figure 5(a) shows the relation of the scale widths between the downstream and upstream sides for each filament. Although there is a large scatter, $`w_\mathrm{u}`$ is systematically smaller than $`w_\mathrm{d}`$ in both the non-thermal and thermal emissions.

Figure 5(b) shows the relation between $`w^\mathrm{s}`$ and $`w^\mathrm{h}`$. We find that $`w_\mathrm{d}^\mathrm{s}`$ is significantly larger than $`w_\mathrm{d}^\mathrm{h}`$, whereas $`w_\mathrm{u}^\mathrm{s}`$ and $`w_\mathrm{u}^\mathrm{d}`$ are comparable with each other. The mean values are $`\overline{w_\mathrm{u}^\mathrm{s}}=4.3`$ arcsec = 0.04 pc, $`\overline{w_\mathrm{d}^\mathrm{s}}=1.3\times 10^2`$ arcsec = 1.1 pc, $`\overline{w_\mathrm{u}^\mathrm{h}}=5.1`$ arcsec = 0.04 pc, and $`\overline{w_\mathrm{d}^\mathrm{h}}=18`$ arcsec = 0.2 pc. Note that the minimum value of $`w_\mathrm{u}^\mathrm{h}`$ is only 0.98 arcsec = 0.01 pc (filament 4). Maybe since their wide scatter is due to that these are the projected values of the possible sheet-like with wavy and/or curved shape.

## 4 Discussion

### 4.1 Thermal Plasma

Although the best-fit NEI model in Table 1 is rejected statistically, it globally fits the thermal emission of SN 1006 as shown in Figure 3. The temperature ($`kT`$=0.24 keV) is similar to the results obtained by Vink et al. (2000), Dyer et al. (2001), and Allen, Petre, & Gotthelf (2001), but lower than that by Koyama et al. (1995). Since the spatial resolution of Chandra enables us to remove the non-thermal photons from the thermal emission more accurately than the ASCA case, the present results should give a more precise description on the thermal plasma. Like the previous observations (Koyama et al., 1995; Allen, Petre, & Gotthelf, 2001), heavy elements, in particular iron, are overabundant, which implies that the X-ray emitting thermal plasma is dominated by the ejecta from type Ia SN. The fact that the thermal emission is enhanced at the inner shell region also suggests the ejecta origin.

From the emission measure ($`E.M.`$) and assuming a uniform density plasma of a prorate shape with the 3-axis radii of 140, 120, and 120 arcsec, we estimate the density $`n_\mathrm{e}`$ in the inner region to be 0.36 cm-3. Then from the best-fit ionization parameter, $`\tau `$ is $`2.6\times 10^{10}`$ s = $`8.3\times 10^2`$ yr, roughly consistent with the age of SN 1006.

Even in the soft2 band profiles, the thermal components are not prominent (see Figure 4), which prevents us from high quality study for the morphology of the thermal plasma. Nevertheless, we found that the profiles of thermal filaments are largely anti-symmetric. The scale width $`w_\mathrm{u}^\mathrm{s}`$ is very sharp and comparable to $`w_\mathrm{u}^\mathrm{h}`$, whereas $`w_\mathrm{d}^\mathrm{s}`$ shows a relatively large scale width. Although $`w_\mathrm{u}^\mathrm{s}`$ couples to $`x_0`$, which is frozen to the best-fit value in the hard band (see Figure 5(b)), we can say that the thermal shock has very sharp rise in upstream.

The downstream scale width $`w_\mathrm{d}^\mathrm{s}`$ is comparable to the shock width derived from the Sedov solution of about 75 arcsec = 0.7 pc for SN 1006 having a 15 arcmin = 8 pc radius (Green, 2001). Therefore, the thermal filaments may nicely trace the density profiles of the Sedov solution.

### 4.2 Non-Thermal Filaments

In this section, we interpret the scale width of hard band X-rays in upstream, $`w_\mathrm{u}^\mathrm{h}`$.

The hardest spectral regions, the non-thermal filaments, are the most probable site of the maximum energy acceleration. From the observation alone, we can not judge whether the filaments are strings or sheets in edge-on configuration. Hester (1987) suggested that thin sheet-like shock fronts are seen as filaments on the edge of the SNR. The filaments seen in figures of Hester (1987) resemble the structure in SN 1006. Furthermore, filaments should be also seen in the inner part of the shell if they are strings, however it is not the observed case. Therefore, we assume that the filaments have sheet-like structure normal to the shock direction. The depth of the sheet is unclear but would be similar to or smaller than the length of the filament. We thus assume the depth of the sheet to be about 1 pc in the following discussion.

The most likely scenario of cosmic ray acceleration at the SNR shock is diffusive shock acceleration (DSA). The predicted results from this model however are highly dependent on many parameters of the magnetic field, such as the angle between the magnetic field and shock normal direction, magnetic field strength and its fluctuation. Here, we investigate the observed profiles based on a DSA model, and estimate the physical quantities such as diffusion coefficient, the direction of magnetic field, maximum energy of accelerated electrons, and the injection efficiency. Although a realistic condition may be more complex as is discussed by Ellison, Berezhko, & Baring (2000) and Berezhko et al. (2002), we at first apply a simplest DSA model that the back reaction of accelerated particles is neglected.

#### 4.2.1 Diffusion coefficient in upstream

The scale widths in the upstream side are largely scattered in the range of 0.98–11 arcsec (Table 2). Since this large scatter would be due to the projected effect of possible sheet-like structure with wavy and/or curved shape, real widths should be smaller than the observed (projected) width. To be conservative, we however adopted the mean value of 4.3 arcsec, or 0.04 pc at 1.8 kpc distance in the following discussion.

We use the results of srcut, the fit of the wide band spectra from X-ray to radio band. The high energy roll-off $`\nu _{\mathrm{rolloff}}`$ srcut is the consequent of either age-limit (acceleration dominant), synchrotron loss, or diffusive escape (see § A.1). The best-fit $`\nu _{\mathrm{rolloff}}`$ at the filaments is $`2.6_{0.7}^{+0.7}\times 10^{17}`$ Hz, which can be converted to $`1.1_{0.3}^{+0.3}`$ keV. These are consistent with the ASCA result of 3.0$`{}_{0.2}{}^{+0.1}\times 10^{17}`$ Hz or $`1.24_{0.08}^{+0.04}`$ keV by Dyer et al. (2001).

Since most of the non-thermal X-ray photons are observed in the downstream region, the synchrotron radiation is mainly due to the downstream region. Using the best-fit $`\nu _{\mathrm{rolloff}}`$, we constrain the maximum energy of electrons $`E_{\mathrm{max}}`$ and magnetic field in downstream $`B_\mathrm{d}`$ from the eq.(A6);

$$E_{\mathrm{max}}B_\mathrm{d}^{0.5}=0.37_{0.06}^{+0.04}[\mathrm{ergs}\mathrm{G}^{0.5}].$$
(2)
In the case of the strong shock, the magnetic field in downstream $`B_\mathrm{d}`$ can be related to that in upstream $`B_\mathrm{u}`$ as

$$B_\mathrm{d}=\left(\mathrm{cos}^2\theta _\mathrm{u}+r^2\mathrm{sin}^2\theta _\mathrm{u}\right)^{\frac{1}{2}}B_\mathrm{u},$$
(3)
where $`r`$ and $`\theta _\mathrm{u}`$ are the compression ratio and the magnetic field angle to the shock normal direction in upstream.

The diffusion coefficients in upstream ($`K_\mathrm{u}`$) is estimated from eq.(A9) ($`u_\mathrm{u}=u_\mathrm{s}`$) as following;

$`K_\mathrm{u}`$ $``$ $`w_\mathrm{u}^\mathrm{h}u_\mathrm{u}=\overline{w_\mathrm{u}^\mathrm{h}}u_\mathrm{s}`$ (4)
$`=`$ $`3.1\times 10^{25}[\mathrm{cm}^2\mathrm{s}^1],`$
where the shock speed $`u_\mathrm{s}`$ is assumed to be 2600 km s-1 (Laming et al., 1996). The diffusion coefficient in upstream is given from Using eq.(A10), we can derive $`K_\mathrm{u}`$ for the electrons of $`E_{\mathrm{max}}`$ as,

$$K_\mathrm{u}=\frac{1}{3}\xi _\mathrm{u}\left(\mathrm{cos}^2\theta _\mathrm{u}+\frac{\mathrm{sin}^2\theta _\mathrm{u}}{1+\xi _\mathrm{u}^2}\right)\frac{E_{\mathrm{max}}}{\mathrm{e}B_\mathrm{u}}c,$$
(5)
where $`\xi _\mathrm{u}(>1)`$ is the fluctuation of the magnetic field.

We can derive following equations from eq.(2), (3), (4), and (5), the maximum energy ($`E_{\mathrm{max}}`$) and magnetic field in upstream ($`B_\mathrm{u}`$) are,

$`E_{\mathrm{max}}`$ $`=`$ $`37\xi _\mathrm{u}^{\frac{1}{3}}\left(\mathrm{cos}^2\theta _\mathrm{u}+{\displaystyle \frac{\mathrm{sin}^2\theta _\mathrm{u}}{1+\xi _\mathrm{u}^2}}\right)^{\frac{1}{3}}`$ (6)
$`\times \left(\mathrm{cos}^2\theta _\mathrm{u}+r^2\mathrm{sin}^2\theta _\mathrm{u}\right)^{\frac{1}{6}}[\mathrm{TeV}],`$
$`B_\mathrm{u}`$ $`=`$ $`4.0\times 10^5\xi _\mathrm{u}^{\frac{2}{3}}\left(\mathrm{cos}^2\theta _\mathrm{u}+{\displaystyle \frac{\mathrm{sin}^2\theta _\mathrm{u}}{1+\xi _\mathrm{u}^2}}\right)^{\frac{2}{3}}`$ (7)
$`\times \left(\mathrm{cos}^2\theta _\mathrm{u}+r^2\mathrm{sin}^2\theta _\mathrm{u}\right)^{\frac{1}{6}}[\mathrm{G}].`$
Since SN 1006 is a type Ia SN located at a high latitude of 460 pc above the Galactic plane ($`l=14.6`$ at 1.8 kpc distance), the interstellar magnetic field would be smaller than the typical value in the Galactic plane of the order of 10 $`\mu `$G. We therefore conservatively assume that $`B_\mathrm{u}`$ = 10 $`\mu `$G, which is slightly larger than the combined results of the X-ray and TeV Gamma-rays (Tanimori et al., 2001). Then eq.(7) becomes,

$$\xi _\mathrm{u}^{\frac{2}{3}}\left(\mathrm{cos}^2\theta _\mathrm{u}+\frac{\mathrm{sin}^2\theta _\mathrm{u}}{1+\xi _\mathrm{u}^2}\right)^{\frac{2}{3}}\left(\mathrm{cos}^2\theta _\mathrm{u}+r^2\mathrm{sin}^2\theta _\mathrm{u}\right)^{\frac{1}{6}}=0.25.$$
(8)
For the parallel magnetic field ($`\theta _\mathrm{u}=0^{}`$), $`\xi _\mathrm{u}`$ becomes smaller than 1. This unrealistically small value of $`\xi _\mathrm{u}`$ conflicts to a usual DSA which assume a parallel magnetic field configuration.

An alternative scenario is oblique magnetic field to the shock normal in this region ($`\theta _\mathrm{u}>0^{}`$). For simplicity, we neglect the back reaction of accelerated particle, then the compression ratio becomes 4 for strong shock limit. The allowed region of $`\theta _\mathrm{u}`$ is given from the eq.(8) as;

$$\theta _\mathrm{u}82^{}.$$
(9)
Thus the magnetic field in upstream near the filaments is almost perpendicular.

Assuming $`\theta _\mathrm{u}=90^{}`$, $`r=4`$, and $`B_\mathrm{u}`$ = 10 $`\mu `$G, $`E_{\mathrm{max}}`$ is given from eq.(6),

$`E_{\mathrm{max}}`$ $`=`$ $`r^{\frac{1}{2}}B_\mathrm{u}^{\frac{1}{2}}0.23_{0.04}^{+0.03}[\mathrm{TeV}]`$ (10)
$`=`$ $`37_7^{+4}[\mathrm{TeV}].`$
When particles are accelerated very efficiently, their back reactions to the shock can not be ignored hence $`r`$ becomes larger than 4. Even if the compression ratio is the largest of $`r=7`$ (Ellison, Berezhko, & Baring, 2000; Berezhko et al., 2002), the allowed range is $`\theta _\mathrm{u}>80^{}`$. Therefore, we can safely predict that the magnetic field in the NE shell of SN 1006 is nearly perpendicular to the shock normal.

As for the downstream region, the observed spatial profile seems to be incompatible with the solution derived by Blandford & Ostriker (1978). The maximum electron energy $`E_{\mathrm{max}}`$ would be determined by the balance of the time scales between the accelerating and the synchrotron cooling. These time scales may depend on the structure and the fluctuation of the magnetic field along the shock normal; both time scales become smaller with larger magnetic field, while the former becomes also smaller with larger fluctuation of the magnetic field.

The shock flow may compresses and partly stretches the magnetic field in the radial direction, which may produce highly disordered magnetic field with small fraction of radial component, as discussed by Reynolds & Gilmore (1993) with the radio polarization data, who reported that only 15–20% of the magnetic-field energy in SN 1006 NE shell in radial polarization, and most of the magnetic field is disordered with the scale smaller than 0.2 pc.

To determine the magnetic field in downstream, we thus need require more complicated processes such as the history of the shock propagation and the non-linear effect, and many assumptions, which is beyond this paper and would leave for a future study.

#### 4.2.2 Injection Efficiency

We define the injection efficiency $`\eta \frac{n_\mathrm{e}^{\mathrm{NT}}}{n_\mathrm{e}^\mathrm{T}}`$, where $`n_\mathrm{e}^{\mathrm{NT}}`$ and $`n_\mathrm{e}^\mathrm{T}`$ are the number densities of non-thermal and thermal electrons in the filaments. The depth of the filament (sheet-like) is assumed to be 1 pc (see § 4.2.1), with uniform electron density.

The non-thermal electron flux (energy) are estimated using the method given in Appendix (§ A.1), where the non-thermal X-ray flux and spectra of each filament are taken from Table 2 and Table 3. We adopt the minimum energy of non-thermal electrons (injection energy) $`E_{\mathrm{min}}`$ to be 0.24 keV (the temperature of the thermal plasma; see Table 1). The magnetic field and the maximum energy are unknown parameters, however we can suggest that the magnetic field in downstream is significantly larger than that in upstream. In this section, we adopt the magnetic field $`B`$ to be perpendicular to the shock normal, $`B_\mathrm{d}=4B_\mathrm{u}`$ = 40 $`\mu `$G arbitrary and the maximum energy $`E_{\mathrm{max}}`$ to be 37 TeV from the discussion in § 4.2.1. The derived number density ($`n_\mathrm{e}^{\mathrm{NT}}`$), the total number ($`N_\mathrm{e}^{\mathrm{NT}}`$), and the total energy ($`E_\mathrm{e}^{\mathrm{NT}}`$) of non-thermal electrons in each filament are summarized in Table 4.

For the estimation of thermal electron flux (energy), we adopted the projected profile of the thermal X-rays (the soft2 band) given in Table 2 and the spectral parameters of thermal plasma given in Table 3. The resultant number density ($`n_\mathrm{e}^\mathrm{T}`$), total number ($`N_\mathrm{e}^\mathrm{T}`$), and total energy ($`E_\mathrm{e}^\mathrm{T}`$) of thermal electrons in each filament are given in Table 4. The injection efficiency $`\eta (\frac{n_\mathrm{e}^{\mathrm{NT}}}{n_\mathrm{e}^\mathrm{T}})`$ is then obtained as is shown in Table 4. All the values $`\eta `$ are nearly identical in each filament of $`1\times 10^3`$ but are about 2 times larger than $`5\times 10^4`$ derived from the ASCA data by Allen, Petre, & Gotthelf (2001), although they also assumed as $`B_\mathrm{d}=40\mu `$G arbitrary. For comparison, we estimate $`\eta `$ with larger magnetic field of $`B_\mathrm{d}=40\mu `$G, and find a larger $`\eta `$ than Allen, Petre, & Gotthelf (2001). They estimated the number density of non-thermal electrons from larger regions than the filaments, because ASCA could not resolve the filaments. On the other hand, we found that the non-thermal electrons are confined in the thin filaments, which are the sites of ongoing acceleration of the non-thermal electrons. We thus regard that the present Chandra result suggests that the injection occurs more locally and more efficiently and as a result, it must be more realistic estimation of $`\eta `$ than that by Allen, Petre, & Gotthelf (2001).

The energy densities of the magnetic field, the thermal plasma, and the non-thermal electrons in the filaments are $`6.4\times 10^{11}`$ ergs cm-3, $`2.6\times 10^{10}`$ ergs cm-3, and $`6.9\times 10^{11}`$ ergs cm-3, respectively. Thus, at the shocked region, the magnetic field and non-thermal electrons are in energy equi-partition and are slightly smaller than the thermal energy (about 30%).

As suggested by Ellison, Berezhko, & Baring (2000), the non-thermal protons should carry larger energy than electrons, hence particle energy becomes larger than that of magnetic and possibly that of thermal. We therefore must consider non-linear effects suggested by Ellison, Berezhko, & Baring (2000) and by Berezhko et al. (2002), however more quantitative scenario including the non-linear effects must be a future work.

## 5 Summary

(1) X-ray emissions from the NE shell of SN 1006 are found to be composites of filaments, clumps and more extended diffuse emissions. We resolved the thermal and non-thermal X-rays using the different spectral shape and morphology.

(2) The spectrum of the thermal component has $`kT`$ = 0.24 keV temperature and relatively small ionization time scale $`n_\mathrm{e}\tau =1.1\times 10^{10}`$ s cm-3. The chemical compositions are overabundant, especially in iron, which suggests that the X-ray emitting plasma originates from ejecta and the progenitor is type Ia.

(3) The non-thermal components can be described with a power-law function with the photon index 2.1–2.3 in the narrow filaments and 2.5 at the inner region of the shell.

(4) The structure of the filaments shows different characteristics in thermal and non-thermal X-rays. The thermal plasma has the scale width of 1 pc in downstream, similar to the shock width derived from Sedov equations. The non-thermal filaments show extremely small scale width of $`w_\mathrm{u}^\mathrm{h}`$ 0.04 pc and $`w_\mathrm{d}^\mathrm{h}`$ 0.2 pc, in upstream and downstream, respectively.

(5) In a diffusive shock acceleration model, the observed thin filaments requires a nearly perpendicular magnetic field with the angle between the magnetic field in upstream (assumed as 10 $`\mu `$G) and the shock normal of larger than $`80^{}`$. The maximum energy of the electrons are 30–40 TeV.

(6) The injection efficiency is estimated to be $`\eta `$ = 1$`\times 10^3`$, suggesting that thermal particles are injected locally and very effectively. Then the energy density of non-thermal electrons becomes comparable to that of the magnetic field and about 30% of the thermal energy density. Thus non-linear effect of the shock structure and acceleration mechanism must be considered.

We thank the anonymous referee for his/her helpful comments. Our particular thanks are due to M. Hoshino, T. Terasawa, T. Yoshida, and S. Inutsuka for their fruitful discussions and comments. A.B. and M.U. are supported by JSPS Research Fellowship for Young Scientists.

## Appendix A Appendix

In this appendix, we briefly introduce the relevant software tools and equations, which are used for the discussions in the text.

### A.1 Electron Spectra and srcut Model

The spectrum of non-thermal electrons accelerated by the diffusive shock is (Bell, 1978);

$`{\displaystyle \frac{\mathrm{d}n_\mathrm{e}^{\mathrm{NT}}}{\mathrm{d}E^{\mathrm{NT}}}}`$ $`=`$ $`\kappa (E^{\mathrm{NT}}+m_\mathrm{e}c^2)(E_{}^{\mathrm{NT}}{}_{}{}^{2}+2m_\mathrm{e}c^2E^{\mathrm{NT}})^{\frac{(p+1)}{2}}\times \mathrm{exp}(E^{\mathrm{NT}}/E_{\mathrm{max}}),`$ (A1)
where $`n_\mathrm{e}^{\mathrm{NT}}`$ and $`E^{\mathrm{NT}}`$ are the number density and the energy of non-thermal electrons, respectively.

The synchrotron radiation power per unit frequency from a single electron of energy $`E`$ in a magnetic field $`B`$ is

$$P(\nu ,\alpha )=\frac{\sqrt{3}\mathrm{e}^3B\mathrm{sin}\alpha }{m_\mathrm{e}c^2}F\left(\frac{\nu }{\nu _c}\right),$$
(A2)
where $`\alpha `$ and $`F(x)`$ are the pitch angle and the function given by Rybicki & Lightman (1979). The peak frequency $`\nu _c`$ is

$$\nu _c=\frac{3cE^2\mathrm{e}B\mathrm{sin}\alpha }{4\pi (m_\mathrm{e}c^2)^3}.$$
(A3)
Convolving the eq.(A1) with (A2), we obtain the spectrum of synchrotron radiation in the pitch angle $`\alpha `$ as,

$$f(\nu ,\alpha )=_{E_{\mathrm{min}}}^{\mathrm{}}dE^{\mathrm{NT}}\frac{\mathrm{d}n_\mathrm{e}^{\mathrm{NT}}}{\mathrm{d}E_\mathrm{e}^{\mathrm{NT}}}P(\nu ,\alpha ).$$
(A4)
Averaged over the pitch angle, we finally obtain the synchrotron radiation energy per unit volume, frequency, and time,

$$f(\nu )=\frac{1}{2}_0^\pi d\alpha \mathrm{sin}\alpha f(\nu ,\alpha ).$$
(A5)
The observed spectrum is fitted with this model spectrum using the Chandra software srcut with 3 free parameters (Reynolds, 1998; Reynolds & Keohane, 1999). The normalization constant and spectral index are determined so as to reproduce the flux and slope in the radio band. These parameters are converted to $`\kappa `$ and $`p`$ in eq.(A1). The other fitting parameter $`\nu _{\mathrm{rolloff}}`$ is a function of $`E_{\mathrm{max}}`$ in eq.(A1) and magnetic field ($`B`$) as is shown below (Reynolds & Keohane, 1999);

$$\nu _{\mathrm{rolloff}}=0.5\times 10^{16}\left(\frac{B}{10\mu \mathrm{G}}\right)\left(\frac{E_{\mathrm{max}}}{10\mathrm{TeV}}\right)^2[\mathrm{Hz}].$$
(A6)
This equation gives constraint on the maximum electron energy $`E_{\mathrm{max}}`$ and magnetic field $`B`$.

Since available radio data of SN 1006 have spatial resolution far larger than the scale of the X-ray filaments, we have no accurate radio flux nor index from the position of the X-ray filaments. We therefore fixed the radio index to the poor resolution radio result of $`p=2.14`$ by Allen, Petre, & Gotthelf (2001) and the radio flux is treated as a free parameter for the present srcut fitting. Allowing the radio index to vary from $`p=2.0`$ to $`p=2.2`$, we see no significant difference in the best-fit parameters within the statistical errors.

### A.2 Diffusion versus Advection

The spatial structure of the relativistic electrons produced by the diffusive shock acceleration across the shock is determined by the competing process of diffusion and advection. The advection time scale ($`\tau _{\mathrm{ad}}`$) is given by

$$\tau _{\mathrm{ad}}=\frac{w}{u},$$
(A7)
where $`w`$ is the scale width of the spatial distribution of the relativistic electrons and $`u`$ is the flow speed.

The diffusion time scale ($`\tau _{\mathrm{dif}}`$) is given from the random-walk theory as;

$$\tau _{\mathrm{dif}}=\frac{w^2}{K},$$
(A8)
where $`K`$ is the diffusion coefficient.

In order that the high energy electrons are accelerated at shock front, electrons in upstream should be diffused back to the shock front against the advection to the downstream side (shock flow). Therefore, $`\tau _{\mathrm{ad}}`$ should be nearly equal to $`\tau _{\mathrm{dif}}`$, hence $`\frac{w_\mathrm{u}}{u_\mathrm{u}}\frac{w_\mathrm{u}^2}{K_\mathrm{u}}`$. We thus obtain;

$$w_\mathrm{u}\frac{K_\mathrm{u}}{u_\mathrm{u}}.$$
(A9)
For more exact formalisms, see e.g. Bell (1978), Blandford & Ostriker (1978), and Drury (1983).

Let the mean free path parallel to the magnetic field be a constant factor $`\xi `$ times the gyro radius $`r_\mathrm{g}`$. Then, the effective diffusion coefficient along the shock normal is given as (Jokipii, 1987; Skilling, 1975)

$`K`$ $`=`$ $`{\displaystyle \frac{1}{3}}\xi r_\mathrm{g}c\left(\mathrm{cos}^2\theta _\mathrm{u}+{\displaystyle \frac{\mathrm{sin}^2\theta _\mathrm{u}}{1+\xi ^2}}\right),`$ (A10)
$`r_\mathrm{g}`$ $`=`$ $`{\displaystyle \frac{E_\mathrm{e}}{\mathrm{e}B}},`$ (A11)
where $`e`$, $`E_\mathrm{e}`$, $`B`$, and $`\theta _\mathrm{u}`$ are the electric charge, the energy of relativistic electrons, and the magnetic field, and the angle between the magnetic field in upstream and the shock normal, respectively.


text:
# Ordered Quantization and the Ehrenfest Time Scale

## I Introduction

Since the earliest days of the quantum theory, the investigation of the differences between the probabilistic concepts in classical Liouville and quantum dynamics has been an important issue. There has been many studies in the last two decades concerning the semi-classical regimen of systems whose classical counterparts exhibits chaos 1 . The question of estimating how long the classical and quantum evolution stay close has been one of the main questions in semi-classical analysis. For classically chaotic flows, there was a conjecture zaslavsky81  which is now rigorously proved hagerdorn00 , that the break time or Ehrenfest time ($`t_E`$) diverges logarithmically with $`\mathrm{}`$. In the classically regular flow, it was suggested in lai  that the behavior of $`t_E`$ with $`\mathrm{}`$ is algebraic ($`\mathrm{}^\delta `$ ), but no universal nature of such behavior has been shown yet.

According to the famous Ehrenfest’s work ehrenfest , for quantum states which are localized enough, the time variation of the mean quantum momentum must be equal to the local force. This statement is exact for quadratic Hamiltonians, but its validity is restricted to a short time scale, the Ehrenfest scale, for more general non-linear systems. Mathematically, the Ehrenfest theorem allows us to write $`𝒪(\widehat{X})=𝒪(\widehat{X})`$ for times smaller than the Ehrenfest time. In this situation, the initial dynamics is described essentially by a mean field approximation, where we have a localized packet obeying classical equations of motion. Then, it is reasonable to expect a rather good agreement between quantum and classical Liouvillian centroids and the classical trajectory. In fact, this scenario have already been reported in the literature ballentine98 .

In this work, we propose a simple analytical scheme to calculate the Ehrenfest time for integrable systems. Our starting point is to propose a classical Liouvillian operator that makes explicit the symmetric form of the usual quantization rules. Using such an operator, we are able to take a general classical function $`𝒪`$ which expresses the time evolution of a physical quantity, quantize it, and evaluate its expectation value in coherent states at each time. The final result, which is shown to be analytically identical to the statistical average calculated through the classical Liouville formalism, is written in terms of a certain differential operator acting onto the original classical function. Since this is just a compact form of expressing the action of the corresponding power series in $`\mathrm{}`$, our recipe automatically leads to a semi-classical expansion around the time evolved classical function. Notice that the equality between quantum and statistical centroids guarantees that we are working within a classical time scale. This allows us to define mathematically the Ehrenfest scale using just the Liouvillian Gaussian average without solving the quantum dynamical problem.

As an example, we present an explicit calculation of the Ehrenfest time for $`N`$-coupled non-linear oscillators, with nonlinearity of order $`2k`$, an even integer number. This model to which we will address the quantum-classical departure issue can be associated to several nonlinearly interacting fields via Kerr-type milburn ; kitagawa  and cross-Kerr type interaction perinova  known to be relevant in quantum optics scully , and also for two quantized vibrational modes of a single trapped ion steinbach . It is an integrable model where the role of nonlinearity can be studied analytically for several quantities, and, in particular, for $`N=k=2`$ the quantum-classical break time has been determined sanz  based on physical properties of the exact quantum states. Although the emphasis of studying the break times in the literature has been on ‘chaotic states’ ballentine ; karkuszewski , here we are concerned with integrable cases, where we are able to derive an analytical expression for the Ehrenfest time.

## II Ordered Quantization

### II.1 Definition

Let us start by presenting a convenient quantization scheme for a single degree of freedom system, which will be easily generalized to degree $`N`$. Consider two non-commuting operators $`\widehat{A}`$ and $`\widehat{B}`$ whose commutator is a $`c`$-number, denoted by $`c`$ ($`[\widehat{A},\widehat{B}]=c`$). Then, a given classical monomial $`a^nb^m`$, where $`a`$ and $`b`$ are canonically conjugated classical variables 111we are considering the cases where $`a`$ and $`b`$ are associated to the operators $`\widehat{A}`$ and $`\widehat{B}`$ trough the coherent state mean values, e.g., $`\alpha _o|\widehat{A}|\alpha _o`$. The canonical pairs $`(\widehat{a},\widehat{a}^{})`$ and $`(\widehat{Q},\widehat{P})`$ are both contempleted in our notation., will be quantized in a symmetric and ordered form, through the prescriptions:

$`a^nb^m𝒮_{\widehat{A},\widehat{B}}\left(\widehat{A}^n\widehat{B}^m\right),`$ (1)
where the super-operator $`𝒮_{\widehat{A},\widehat{B}}`$ is given by

$`𝒮_{\widehat{A},\widehat{B}}=e^{\frac{1}{2}[\widehat{A},\widehat{B}]_{\widehat{A}}_{\widehat{B}}}={\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\left(\frac{1}{2}[\widehat{A},\widehat{B}]\right)^k}{k!}}_{\widehat{A}}^k_{\widehat{B}}^k.`$ (2)
The above index $`(\widehat{A},\widehat{B})`$ corresponds to the ordering with $`\widehat{A}`$ on the left and $`\widehat{B}`$ on the right. Since $`[_{\widehat{A}},_{\widehat{B}}]=[_{\widehat{A}},\widehat{B}]=[\widehat{A},_{\widehat{B}}]=0`$, $`𝒮_{\widehat{A},\widehat{B}}`$ and $`𝒮_{\widehat{B},\widehat{A}}`$ are classical unitary differential super-operators, satisfying $`𝒮_{\widehat{A},\widehat{B}}𝒮_{\widehat{B},\widehat{A}}=𝒮_{\widehat{B},\widehat{A}}𝒮_{\widehat{A},\widehat{B}}=1`$. The prescription defined in expression (1) leads to different orderings (depending on which variable is chosen to be $`a`$ or $`b`$) for the same original classical function and, therefore, the associated symmetric operators must be the same. From this consideration, we deduce the following ordering formulas:

$$\begin{array}{c}\widehat{A}^n\widehat{B}^m=𝒮_{\widehat{B},\widehat{A}}^2\widehat{B}^m\widehat{A}^n=e^{+[\widehat{A},\widehat{B}]_{\widehat{A}}_{\widehat{B}}}\widehat{B}^m\widehat{A}^n,\hfill \\ \\ \widehat{B}^m\widehat{A}^n=𝒮_{\widehat{A},\widehat{B}}^2\widehat{A}^n\widehat{B}^m=e^{[\widehat{A},\widehat{B}]_{\widehat{A}}_{\widehat{B}}}\widehat{A}^n\widehat{B}^m.\hfill \end{array}$$
(3)
Using a classical displacement operator, $`e^{a_x}f(x)=f(x+a)`$, which can be used to write $`e^{\frac{c}{2}_{\widehat{A}}_{\widehat{B}}}\widehat{A}^n\widehat{B}^m=(\widehat{A}\frac{c}{2}_{\widehat{B}})^n\widehat{B}^m`$, one can show that these results reproduce those in the textbooks (see e.g. Louisell louisell ). A simple example shows that our quantization scheme leads to the usual quantization rules, but in an automatically ordered form. Consider the following product of classical canonical phase space variables: $`q^2p`$. According to the usual rules, this is transformed in a totally symmetric operator cahill  of the form $`\frac{1}{3}(\widehat{Q}^2\widehat{P}+\widehat{P}\widehat{Q}{}_{}{}^{2}+\widehat{Q}\widehat{P}\widehat{Q})`$; and using the commutation relation $`[\widehat{Q},\widehat{P}]=ı\mathrm{}`$, re-written as $`\widehat{Q}{}_{}{}^{2}\widehat{P}ı\mathrm{}\widehat{Q}`$. But this is exactly the result produced by the expression in (1) with the choice $`(a,b)=(q,p)`$, just by making some few derivatives. We finally note that the super-operators $`𝒮_{\widehat{A},\widehat{B}}`$ has already appeared in the literature in the case of the canonical phase space operators $`(\widehat{A},\widehat{B})=(\widehat{Q},\widehat{P})`$, in slightly different contexts royer91 .

### II.2 Normal Ordering in Bosonic Operators

As an immediate application of the formulas presented above, we will express the totally symmetric ordered expression $`𝒮_{\widehat{Q},\widehat{P}}\widehat{Q}^n\widehat{P}^m`$ in a form more suitable for our purpose of taking expectation values in the Weyl-Heisenberg coherent states. We start by expressing such a polynomial in terms of bosonic creation and annihilation operators

$`𝒮_{\widehat{Q},\widehat{P}}\widehat{Q}^n\widehat{P}^m`$ $`=`$ $`e^{\frac{ı\mathrm{}}{2}_{\widehat{Q}}_{\widehat{P}}}\widehat{Q}^n\widehat{P}^m`$
$`=`$ $`z_{n,m}e^{\frac{1}{4}\left(_{\widehat{a}}^2_{\widehat{a}^{}}^2\right)}\left(\widehat{a}+\widehat{a}^{}\right)^n\left(\widehat{a}\widehat{a}^{}\right)^m,`$
where $`z_{nm}=(ı)^m\left(\sqrt{\mathrm{}/2}\right)^{n+m}`$. Now, using the fact that

$`\left(\widehat{a}\pm \widehat{a}^{}\right)^n`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{n}{}}}\left(\pm 1\right)^{nk}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)e^{\frac{1}{2}_{\widehat{a}}_{\widehat{a}^{}}}\left(\widehat{a}^{nk}\widehat{a}^k\right),`$ (5)
we can re-express (II.2) as

$`𝒮_{\widehat{Q},\widehat{P}}\widehat{Q}^n\widehat{P}^m={\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \underset{l=0}{\overset{m}{}}}(1)^{ml}z_{n,m}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{m}{l}}\right)\times `$
$`\times e^{\frac{1}{4}\left(_{\widehat{a}}^2_{\widehat{a}^{}}^2\right)}e^{\frac{1}{2}_{\widehat{a}_2}_{\widehat{a}_1^{}}}e^{\frac{1}{2}_{\widehat{a}_4}_{\widehat{a}_3^{}}}\left(\widehat{a}_1^{nk}\widehat{a}_2^k\widehat{a}_3^{ml}\widehat{a}_4^l\right),`$ (6)
with the sub-indexes that we introduced, to indicate where the action of the differentiation should take place. At the end of calculation we must erase all these sub-indexes. Using the relation (3), we get

$`𝒮_{\widehat{Q},\widehat{P}}\widehat{Q}^n\widehat{P}^m={\displaystyle \underset{k=0}{\overset{n}{}}}{\displaystyle \underset{l=0}{\overset{m}{}}}(1)^{ml}z_{n,m}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{m}{l}}\right)\times `$ (7)
$`\times e^{\frac{1}{4}\left(_{\widehat{a}}^2_{\widehat{a}^{}}^2\right)}e^{\frac{1}{2}_{\widehat{a}_2}_{\widehat{a}_1^{}}}e^{\frac{1}{2}_{\widehat{a}_4}_{\widehat{a}_3^{}}}e^{_{\widehat{a}_2}_{\widehat{a}_3^{}}}\left(\widehat{a}_1^{nk}\widehat{a}_3^{ml}\widehat{a}_2^l\widehat{a}_4^k\right).`$
This is the normal ordered expression (in the creation and annihilation operators) for the original $`\widehat{Q}\widehat{P}`$-ordered monomial. Now, it is a simple matter to calculate its expectation value.

### II.3 Coherent States Representation

Since we are interested in the connection between quantum and classical mechanics, the coherent state basis appears as the most appropriated one. In particular, it will be of particular interest for us to evaluate the expectation value in coherent states of some operator products like the one treated in previous section. Then, we first calculate the matrix elements in the coherent state basis of the operator function given in Eq.(7)

$`{\displaystyle \frac{\alpha _1|𝒮_{\widehat{Q},\widehat{P}}\widehat{Q}^m\widehat{P}^m|\alpha _2}{\alpha _1|\alpha _2}}=z_{n,m}e^{\frac{1}{4}\left(_{\alpha _2}^2_{\alpha _1^{}}^2\right)}\{e^{_a_d}e^{\frac{1}{2}_a_b}\times `$
$`\times e^{\frac{1}{2}_c_d}\left[\sqrt{{\displaystyle \frac{\mathrm{}}{2}}}(a+b)\right]^n\left[\sqrt{{\displaystyle \frac{\mathrm{}}{2}}}{\displaystyle \frac{(cd)}{ı}}\right]^m\}_{\genfrac{}{}{0pt}{}{a=c=\alpha _2}{b=d=\alpha _1^{}}}.`$ (8)
Setting now $`\alpha _1=\alpha _2=\alpha _0=\frac{q_0+ıp_0}{\sqrt{2\mathrm{}}}`$ and performing adequate variable transformations, we finally get

$`\alpha _0|𝒮_{\widehat{Q},\widehat{P}}\widehat{Q}^n\widehat{P}^m|\alpha _0=e^{\frac{\mathrm{}}{4}_0^2}q_0^np_0^m,`$ (9)
where

$$\begin{array}{c}_0{}_{}{}^{2}=_{q_0}{}_{}{}^{2}+_{p_0}{}_{}{}^{2},\\ \\ (q_0,p_0)=(\alpha _0|\widehat{Q}|\alpha _0,\alpha _0|\widehat{P}|\alpha _0).\end{array}$$
(10)
The expectation value given in (9) must be calculated through the series expansion of the classical operator $`e^{\frac{\mathrm{}}{4}_0^2}`$, which gives a natural expansion in powers of $`\mathrm{}`$, showing its semi-classical nature. In fact, in the classical limit $`\mathrm{}0`$, the quantum expectation value of the operator function reduces to a purely classical function calculated at the center of the coherent packet located at $`(q_0,p_0)`$.

Another interesting relation can be obtained from a similar calculation

$$\begin{array}{c}\widehat{Q}^n\widehat{P}^m=e^{\frac{\mathrm{}}{4}_0^2}\left(e^{\frac{ı\mathrm{}}{2}_{q_0}_{p_0}}q_0^np_0^m\right),\end{array}$$
(11)
which implies that

$$\begin{array}{c}\frac{1}{ı\mathrm{}}[\widehat{Q}^n,\widehat{P}^m]=e^{\frac{\mathrm{}}{4}_0^2}\left(\frac{\mathrm{sin}\frac{\mathrm{}}{2}_{q_0}_{p_0}}{\frac{\mathrm{}}{2}}\right)q_0^np_0^m.\end{array}$$
(12)
It is important to note that the term within the parentheses in Eq.(11) is exactly the Weyl transform of the operator $`\widehat{Q}^n\widehat{P}^m`$ royer91 . The extra operator factor $`e^{\frac{\mathrm{}}{2}_0^2}`$ is what makes reference to the width of the coherent packet, as we shall see latter. These results also point out to the existence of an asymptotic classical limit ($`\mathrm{}0`$) for such expectation values in coherent states.

## III Short Time Quantization

The usual quantization rules are defined in the Heisenberg picture, where the Hamilton’s equations solution $`q_{}(q,p,t)`$ and $`p_{}(q,p,t)`$, are transformed into Heisenberg operators $`\widehat{Q}_H(\widehat{Q},\widehat{P},t)`$ and $`\widehat{P}_H(\widehat{Q},\widehat{P},t)`$, where $`\widehat{Q}`$ and $`\widehat{P}`$ denote Schrödinger operators. On the other hand, since the Heisenberg and Schrödinger pictures coincide at $`t=0`$, any scheme of quantization based on the Schrödinger picture would describe reasonably the quantum world for very short times. However, as long as we are interested in analyzing the quantum operator evolution only during a classical (mean field) time scale, this would suffice.

In this context, consider a classically time evolved function $`𝒪`$ that has the following well defined expansion

$`𝒪(q,p,t)={\displaystyle \underset{n,m=0}{\overset{\mathrm{}}{}}}c_{n,m}(t)q^np^m.`$ (13)
The coefficients $`c_{n,m}(t)`$ contain all the time dependence. Now, applying the operator of ordered quantization $`𝒮_{\widehat{Q},\widehat{P}}`$ and taking the average in coherent states we obtain

$`\alpha _0|\widehat{𝒪}(\widehat{Q},\widehat{P},t)|\alpha _0=e^{\frac{\mathrm{}}{4}_0^2}𝒪(q_0,p_0,t),`$ (14)
where we have used the result (9). In order to estimate for how long the result (14) could be trusted, we present the calculation of the classical statistical counterpart of the problem. In the classical Liouvillian formalism, the following mean value is defined in the phase space $`(q,p)`$

$$𝒪(q_0,p_0,t)=𝑑q𝑑p\rho (q,p,t)𝒪(q,p,0),$$
(15)
where ($`q_0`$,$`p_0`$) stands for the center of the initial distribution $`\rho (q,p,0)`$. Then, since $`\rho (q,p,t)=e^t\rho (q,p,0)`$ and $`𝒪(q,p,t)=e^t𝒪(q,p,0)`$, taking into account the fact that the volume in the phase space is preserved, we can re-write the previous equation in a Heisenberg-like form

$`𝒪(q_0,p_0,t)={\displaystyle 𝑑q𝑑p\rho (q,p,0)𝒪(q,p,t)}.`$ (16)
Consider now a Gaussian initial distribution with width $`\sigma `$, and centered at the point $`(q_0,p_0)`$. Performing the variable transformations, $`(qq_0)=x`$ and $`(pp_0)=y`$, and using the classical displacement operator again, we can rewrite Eq.(16) as

$`𝒪={\displaystyle 𝑑x𝑑y\frac{e^{\frac{x^2}{\sigma }}}{\sqrt{\pi \sigma }}\frac{e^{\frac{y^2}{\sigma }}}{\sqrt{\pi \sigma }}e^{x_q}e^{y_p}𝒪(q,p,t)}.`$ (17)
Notice that the function $`𝒪`$ inside the integral no longer depends on the integration variables. Hence, by completing squares and performing formally the integration, we finally obtain

$`𝒪(q_0,p_0,t)=e^{\frac{\sigma }{4}_0^2}𝒪(q_0,p_0,t).`$ (18)
Now, we have proved that the effect of applying the operator $`e^{\frac{\sigma }{4}_0^2}`$ is exactly that of smoothing the function $`𝒪`$ through a Gaussian mean, where $`\sigma `$ is related to the width of the Gaussian distribution to be used in the smoothing. Accordingly, this operator has already been used to express the $`Q`$-function as the Gaussian smoothing of the Wigner function cessa . It is important to notice that this result is exact for Gaussian statistical averages, i.e., it was derived from first principles without any approximation. The only implicit assumption used in the calculation was the existence of the derivatives in all orders for the function $`𝒪`$.

The comparison between expressions (14) and (18), with $`\sigma =\mathrm{}`$, confirms the fact that our proposal of a Schrödinger quantization for classical function should be adequate either at a short time scale or in cases in which the classical function depends linearly on the phase space coordinates $`(q_0,p_0)`$ (e.g. harmonic oscillator). The results (14) and (18) can easily be extended for higher degrees of freedom.

At this point one might formulate a question. First, we interpret the action of the operator $`e^{\frac{\mathrm{}}{4}_0^2}`$ on any classical function, as an exact factorization of the effect of a Gaussian wave packet contribution. In fact, for completely localized distribution ($`\mathrm{}=0`$) there are no corrections. Then, one may ask: is it possible, in the particular case of coherent separable initial states, to express the exact quantum centroid, $`\widehat{𝐑}(t)=(\widehat{Q}(t),\widehat{P}(t))`$ in terms of the smoothed form $`e^{\frac{\mathrm{}}{4}^2}𝐫_𝐜(t)`$, where $`𝐫_𝐜(t)=(q_c(t),p_c(t))`$ is a certain classical dynamical vector in phase space? The answer will be positive if we are able to calculate the inverse of the Gaussian smoothing operation on the quantum centroid vector in the phase space

$`𝐫_𝐜(t)=e^{\frac{\mathrm{}}{4}_0^2}\alpha _0|\widehat{𝐑}(t)|\alpha _0.`$ (19)
This function $`𝐫_𝐜(t)`$ will then be a coherent quantum trajectory, in the sense that it will carry all quantum dynamical information possible to put in a trajectory, except the contribution due to a Gaussian smoothing effect. For the trivial case of $`N`$ non-interacting harmonic oscillators and the case of a bilinear coupling between two harmonic oscillators, the quantum coherent and classical trajectories coincide ( $`\mathbf{(}q_c(t),p_c(t)\mathbf{)}=\mathbf{(}q_0(t),p_0(t)\mathbf{)}`$ ), as expected.

## IV The Ehrenfest Time

### IV.1 Formal Definition

Now, we have all the necessary tools to undertake the problem of estimating the Ehrenfest time in the case of an initially separable coherent Gaussian wave packet. Assuming that the break instant occurs when the first order correction in $`\mathrm{}`$ becomes as important as the original vector, we apply the smoothing process to the phase space classical trajectory vector of a system with $`N`$ degrees of freedom $`𝐫(t)=\mathbf{(}q_1(t),p_1(t),\mathrm{},q_N(t),p_N(t)\mathbf{)}`$. Expanding the smoothing operator up to the first order in $`\mathrm{}`$, we obtain

$`\widehat{𝐑}(t)=e^{\frac{\mathrm{}}{4}\mathbf{}^2}𝐫(t)𝐫(t)+{\displaystyle \frac{\mathrm{}}{4}}\mathbf{}^2𝐫(t),`$ (20)
where $`\mathbf{}^2=_{i=1}^N_i^2`$ is the $`N`$-dimensional Laplacian operator. Now, we formally define the Ehrenfest time $`t_E`$ as being the instant at which the magnitude of the difference between the quantum centroid and the corresponding classical vector in phase space becomes equal to the magnitude of the initial classical vector. Mathematically this condition can be expressed as follows:

$`{\displaystyle \frac{\widehat{𝐑}(t_E)𝐫(t_E)}{𝐫(0)}}=1.`$ (21)
Finally, using Eq.(20), we obtain the prescription for the analytical calculation of the Ehrenfest time

$`{\displaystyle \frac{\mathrm{}}{4}}{\displaystyle \frac{\mathbf{}^2𝐫(t_E)}{𝐫(0)}}=1.`$ (22)
It is remarkable that we do not need to solve the quantum equations of motion to find $`t_E`$. In what follows, we will calculate the above defined break time for the system of $`N`$ coupled non-linear oscillators.

### IV.2 Application to a Non-Linear System

Consider the following classical Hamiltonian

$`={\displaystyle \underset{i=1}{\overset{N}{}}}\omega _i\left({\displaystyle \frac{p_i^2+q_i^2}{2}}\right)+g\left[{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{p_i^2+q_i^2}{2}}\right)\right]^k,`$ (23)
where $`k1`$ is integer and $`g`$ is the only coupling constant of the system, and from which we define the characteristic classical action $`\mathrm{\Lambda }=\underset{i=1}{\overset{N}{}}\left(\frac{p_i^2+q_i^2}{2}\right)`$ remark . The equations of motion can be solved by noticing that $`\mathrm{\Lambda }`$ itself is a constant of motion. The result reads

$`\left[\begin{array}{c}q_i(t)\\ \\ p_i(t)\end{array}\right]=\left[\begin{array}{cc}\mathrm{cos}\mathrm{\Theta }_i(t)& \mathrm{sin}\mathrm{\Theta }_i(t)\\ & \\ \mathrm{sin}\mathrm{\Theta }_i(t)& \mathrm{cos}\mathrm{\Theta }_i(t)\end{array}\right]\left[\begin{array}{c}q_i\\ \\ p_i\end{array}\right],`$ (33)
where $`(q_i,p_i)`$ are the initial conditions of the $`i`$-th oscillator and $`\mathrm{\Theta }_i(t)=\omega _it+gtk\mathrm{\Lambda }^{k1}`$ is a rotation angle in phase space. Finally, using (22) and (33), we get for the Ehrenfest time of this system

$`t_E`$ $``$ $`{\displaystyle \frac{1}{k(k1)}}\left[{\displaystyle \frac{1}{g\mathrm{\Lambda }^{k1}}}\left({\displaystyle \frac{2\mathrm{\Lambda }}{\mathrm{}}}\right)^{0.5}\right]\left(1{\displaystyle \frac{\mathrm{}k^2}{8\mathrm{\Lambda }}}\right),`$ (34)
where we kept only the first order terms in $`\mathrm{}`$. The cases $`(N=1,k=2)`$ and $`(N=2,k=2)`$ can be seen to reproduce respectively the results calculated in ref.berman  and carolina . Trivial limits are also contemplated by the above result, namely, the cases of an harmonic system ($`g=0`$ or $`k=1`$), for which $`t_E\mathrm{}`$. Moreover, our result is in accordance with some conjectures predicting the general form $`\frac{1}{\mathrm{\Omega }}\left(\frac{S}{\mathrm{}}\right)^\delta `$ for the break time of classically integrable systems zurek , where $`\mathrm{\Omega }=\frac{1}{k(k1)g\mathrm{\Lambda }^{k1}}`$ and $`S2\mathrm{\Lambda }`$ are respectively the typical frequency and classical action of the system in consideration. We illustrate in figure 1 the Ehrenfest scale predicted by the expression (34) for the case $`(N=2,k=2)`$, where we can see an algebraic ($`t^2`$) short time departure. We also note that the first order correction to the Ehrenfest time scale seems to indicate a more appropriated parameter to measure the classicallity in this system: $`\frac{\mathrm{}k^2}{8\mathrm{\Lambda }}1`$.

## V Conclusion

In this paper we proposed a scheme for an ordered symmetric quantization based on the action of certain differential operators. As an application, we quantized a classically time evolved dynamical functions of the canonical variables in phase space, and showed that such a procedure is adequate during the Ehrenfest time scale for separable coherent initial states. This allowed us to propose a formal definition for the Ehrenfest time in terms of the phase space Laplacian operator acting on the classical solutions of the equations of motion. We performed an explicit calculation for a system of $`N`$ coupled nonlinear ($`2k`$-th order) oscillators and calculated the Ehrenfest time for general $`N1`$ and $`k1`$. The results are shown to agree with the results known in the literature for some particular cases.

###### Acknowledgements.

It is a pleasure to acknowledge M.C. Nemes and A.C. Oliveira for many helpful discussions, and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (Contracts No.300651/85-6,, No.146010/99-0) and Fundação de Amparo à Pesquisa de São Paulo (FAPESP) (Contract No.98/13617-4) for financial support.