format: text

text:
# A scale-relativistic derivation of the Dirac Equation

## I Introduction

The theory of scale relativity generalizes to scale transformations the principle of relativity, which has been applied by Einstein to motion laws. It is based on the giving up of the assumption of space-time coordinate differentiability, which is usually retained as an implicit hypothesis in current physics. Even though this hypothesis can be considered as roughly valid in the classical domain, it is clearly broken by the quantum mechanical behavior. It has indeed been pointed out by Feynman (see, e.g., Ref. ) that the typical paths of quantum mechanics are continuous but non-differentiable. In the present paper, we apply the scale-relativistic approach to the understanding of the nature of the bi-spinors and of the Dirac equation in a space-time representation.

One of the fundamental reasons for jumping to a scale-relativistic description of nature is to increase the generality of the first principles retained. Einstein’s principle of relativity relies on the assumption that coordinates are a priori differentiable. It is therefore not fully general, since it applies to coordinate transformations that are continuous and, at least two times, differentiable. The aim of the theory of scale relativity is to bring to light laws and structures that would be the manifestation of more general transformations, namely, continuous ones, either differentiable or not. The standard “general relativistic” theory will be automatically recovered as the special differentiable case, while new effects would be expected from the non-differentiable part.

Giving up the assumption of differentiability has important physical consequences: one can show  that spaces of topological dimension $`D_T`$, which are continuous but non-differentiable, are characterized by a $`D_T`$ measure which becomes explicitly dependent on the resolution (i.e., the observation scale) $`ϵ`$ at which it is considered and tends to infinity when the resolution interval $`ϵ`$ tends to zero. Therefore, a non-differentiable space-time continuum is necessarily fractal, in the general meaning given to this word by Mandelbrot , who thus means, not only self-similar, but explicitly scale-dependent in the more general way. This statement naturally leads to the proposal of a geometric tool adapted to construct a theory based on such premises, namely, a fractal space-time (see, e.g., Refs. ). Such an explicit and fundamental dependence of physical laws on the internal “resolution” scales leads to constrain the new scale laws (that need to be constructed) by a “principle of scale relativity” . This principle extends to scale transformations (i.e., contractions and dilations of resolutions) the Galileo-Einstein principle of relativity, that was up to now applied only to displacements and motion. This allows to include the resolutions in the definition of the state of the reference system and to require scale covariance of the equations of physics under scale transformations. A main consequence is that the geodesics of such a non-differentiable space-time are themselves fractal and in infinite number between any two points. This strongly suggests that the description could borrow some tools from the mathematical formalism of statistics, i.e., will provide a probabilistic interpretation, therefore naturally extending the realm of quantum behaviour to a larger spectrum of available scales. The domains of application are typically the asymptotic scale ranges of physics: (i) Microphysics: small length-scales and time-scales. (ii) Cosmology: large length-scales. (iii) Complex chaotic systems: long time-scales.

In earlier works  the Schrödinger equation and the (free) Klein-Gordon equation have been established, from the first principles of this theory, as geodesic equations on a fractal space/space-time. The third main evolution equation of standard quantum physics, the (free) Dirac equation is here derived, following the same trend of thought.

As in these earlier works, the general method followed here, in order to reach a genuine demonstration of the Dirac equation from first principles, involves the following steps: (i) introduce a scale-space of internal resolutions, aimed at describing the underlying fractal structures of geodesic trajectories (that are identified with the “particles” in this approach); (ii) perform this description by writing scale differential equations coming under the principle of relativity, then solving them; (iii) express the effects induced on the laws of motion (space-time displacement) by the internal fractal structures of the underlying scale-space in terms of a covariant derivative $`{}_{}{}^{}d\text{´}/ds`$; when applied on the classical position variable, this covariant derivative changes the nature of the classical velocity, which becomes complex (Schrödinger and Klein-Gordon equations) , quaternionic (Pauli equation), then bi-quaternionic (Dirac equation) at various depth levels of the theory; (iv) use the principle of relativity in the form of a generalized equivalence/strong covariance principle in order to write the motion equation as a geodesic equation, i.e. as a free motion equation of the form $`{}_{}{}^{}d\text{´}^2x^\mu /ds^2=0`$; (v) expand the covariant derivative and make a change of variable in order to introduce the quantum tool (the wave function is defined as a re-expression of the action), then integrate and obtain the fundamental equations of quantum mechanics.

In Sec. II, we actualize the former works dealing with this issue, proposing a more accurate and profound interpretation of the nature of the transition from the non-differentiable (fractal scales) to the differentiable (classical scales) domain, and carefully justifying the different key choices made at each main step of the reasoning. Then we derive, in Sec. III, the bi-quaternionic form of the Klein-Gordon equation, from which the Dirac equation naturally proceeds, as shown in Sec. IV. Section. V is devoted to the conclusion.

## II The fractal to classical scale transition revisited

### A Scale invariance and Galilean scale relativity

A power law scale dependence, frequently encountered in natural systems, is described geometrically in terms of fractals , and algebrically in terms of the renormalization group . As we show below, such simple scale-invariant laws can be identified with the “Galilean” version of scale-relativistic laws.

Consider a non-differentiable (fractal) curvilinear coordinate $`(x,ϵ)`$, that depends on some space-time variables $`x`$ and on the resolution $`ϵ`$. Such a coordinate generalizes to non-differentiable and fractal space-times the concept of curvilinear coordinates introduced for curved Riemannian space-times in Einstein’s general relativity . Rather than considering only the strict non-differentiable mathematical object $`(x)`$, we are interested in its different versions at various resolutions $`ϵ`$. Such a point of view is particularly well-adapted to applications in physics, as any real measurement is always performed at a finite resolution. In this framework, $`(x)`$ becomes the limit when $`ϵ0`$ of the family of functions $`(x,ϵ)`$. But while $`(x,0)`$ is non-differentiable, $`(x,ϵ)`$, which we call a fractal function , is differentiable for all $`ϵ0`$. The physics of a given process will therefore be completely described if we can know $`(x,ϵ)`$ for all values of $`ϵ`$. $`(x,ϵ)`$, being differentiable when $`ϵ0`$, can be the solution to differential equations involving the derivatives of $``$ with respect to both $`x`$ and $`ϵ`$.

Let us now apply an infinitesimal dilation $`ϵϵ^{}=ϵ(1+d\rho )`$ to the resolution. Being, at this stage, interested in pure scale laws, we omit the $`x`$ dependence in order to simplify the notation and obtain, at first order,

$$(ϵ^{})=(ϵ+ϵd\rho )=(ϵ)+\frac{(ϵ)}{ϵ}ϵd\rho =(1+\stackrel{~}{D}d\rho )(ϵ),$$
(1)
where $`\stackrel{~}{D}`$ is, by definition, the dilation operator. The identification of the two last members of this equation yields

$$\stackrel{~}{D}=ϵ\frac{}{ϵ}=\frac{}{\mathrm{ln}ϵ}.$$
(2)
This well-known form of the infinitesimal dilation operator shows that the “natural” variable for the resolution is $`\mathrm{ln}ϵ`$, and that the expected new differential equations will involve quantities such as $`(x,ϵ)/\mathrm{ln}ϵ`$. An example of equations describing scale dependence is given by the renormalization group equations, as first proposed by Wilson .

We limit the present approach to the consideration of the simplest form that can be exhibited by this class of equations, leaving to future works the study of more complete cases. Let us built such an equation by stating that the variation of $``$ under an infinitesimal scale transformation $`d\mathrm{ln}ϵ`$ depends only on $``$ itself; namely, $``$ determines the whole physical behavior, including the behavior under scale transformations. We thus write

$$\frac{(x,ϵ)}{\mathrm{ln}ϵ}=\beta ().$$
(3)
As we are interested in the simplest form for such an equation, we expand $`\beta ()`$ in powers of $``$. This can always be done since $``$ may be renormalized, dividing it by its largest value, in such a way that the new variable $``$ remains $`1`$ in its variation domain. We obtain, to the first order, the linear equation

$$\frac{(x,ϵ)}{\mathrm{ln}ϵ}=a+b,$$
(4)
of which the solution is

$$(x,ϵ)=_0(x)\left[1+\zeta (x)\left(\frac{\lambda }{ϵ}\right)^b\right],$$
(5)
where $`\lambda ^b\zeta (x)`$ is an integration constant and $`_0=a/b`$.

For $`b<0`$, Eq. (5) gives a fractal scale-invariant behavior at small scales ($`ϵ\lambda `$), with a fractal dimension $`D=D_Tb`$. For a curvilinear coordinate $``$, the topological dimension $`D_T`$ is equal to $`1`$. The scale dimension, $`\delta =DD_T`$, is defined, following Mandelbrot , as

$$\delta =\frac{d\mathrm{ln}}{d\mathrm{ln}(\lambda /ϵ)}.$$
(6)
When $`\delta `$ is constant, one obtains an asymptotic power law resolution dependence

$$(x,ϵ)=_0(x)\left(\frac{\lambda }{ϵ}\right)^\delta .$$
(7)
Let us now check that for such a simple self-similar scaling law the new scale laws to be constructed come indeed under the principle of relativity extended to scale transformations of the resolutions $`ϵ`$, as we have stressed in the introduction. The structure of the corresponding group of scale transformations can be established from the fact that the involved quantities transform, under a scale transformation $`ϵϵ^{}`$, as

$$\mathrm{ln}\frac{(ϵ^{})}{_0}=\mathrm{ln}\frac{(ϵ)}{_0}+\delta (ϵ)\mathrm{ln}\frac{ϵ}{ϵ^{}},$$
(8)
$$\delta (ϵ^{})=\delta (ϵ).$$
(9)
These transformations have exactly the structure of the Galileo group, as confirmed by the dilation composition law, $`ϵϵ^{}ϵ^{\prime \prime }`$, which writes

$$\mathrm{ln}\frac{ϵ^{\prime \prime }}{ϵ}=\mathrm{ln}\frac{ϵ^{}}{ϵ}+\mathrm{ln}\frac{ϵ^{\prime \prime }}{ϵ^{}}.$$
(10)
It is worth noting that Eq. (5) gives also a transition from a fractal to a non-fractal behavior beyond some transition scale $`\lambda `$.

### B Transition from non-differentiability (fractal scales) to differentiability (classical scales)

The first consequence of the giving up of the coordinate differentiability is the differential proper time symmetry breaking. Relativistic motion involves a four-dimensional space-time with the proper time $`s`$ as a curvilinear parameter. Strictly, the non-differentiability of the coordinates implies that the four-velocity

$$V^\mu =\frac{dX^\mu }{ds}=\underset{ds0}{lim}\frac{X^\mu (s+ds)X^\mu (s)}{ds}$$
(11)
is undefined. This means that, when $`ds`$ tends to zero, either the ratio $`dX^\mu /ds`$ tends to infinity, or it fluctuates without reaching any limit. However, as recalled in the introduction, continuity and non-differentiability imply an explicit scale dependence of the various physical quantities, and therefore of the velocity, $`V^\mu `$. We thus replace the differential, $`ds`$, by a scale variable, $`\delta s`$, and consider $`V^\mu (s,\delta s)`$ as an explicit fractal function of this variable. For a constant fractal dimension $`D`$, the resolution in $`X^\mu `$, $`ϵ_\mu (\delta s)`$, corresponds to the resolution in $`s`$, $`\delta s`$, according to

$$ϵ_\mu \delta s^{1/D}.$$
(12)
The advantage of this method is that, for any given value of the resolution, $`\delta s`$, differentiability in $`s`$ is recovered, which allows to use the differential calculus, even when dealing with non-differentiability.

The scale dependence of the velocity suggests that we complete the standard equations of physics by new differential equations of scale. We therefore apply to the velocity and to the differential element, now interpreted as a resolution, the reasoning applied to the fractal function $``$ in Sec. II A. Writing the simplest possible equation for the variation of the velocity $`V^\mu (s,ds)`$ in terms of the new scale variable $`ds`$, as a first order renormalization group-like differential equation of the form of Eq. (3); then, Taylor-expanding it, as in Eq. (4), using the fact that $`V^\mu <1`$ (in motion-relativistic units $`c=1`$), we obtain the solution

$$V^\mu =v^\mu +w^\mu =v^\mu \left[1+\zeta ^\mu \left(\frac{\tau }{ds}\right)^{11/D}\right],$$
(13)
where we have set $`b=1/D1`$. Here, $`v^\mu `$ is the “classical part” of the velocity, $`\overline{𝒞\mathrm{}}V^\mu `$, (see below the definition of the classical part operator $`\overline{𝒞\mathrm{}}`$), $`w^\mu `$ is the explicitly scale-dependent “fractal part” and $`\tau `$ and $`\zeta ^\mu `$ are chosen such that $`\overline{𝒞\mathrm{}}\zeta ^\mu =0`$ and $`\overline{𝒞\mathrm{}}(\zeta ^\mu )^2=1`$. We recognize here the combination of a typical fractal behavior, with a fractal dimension $`D`$, and of a breaking of the scale symmetry at the scale transition $`\tau `$, which is identified with the Compton scale of the system ($`\tau =\mathrm{}/mc`$), since $`V^\mu v^\mu `$, when $`ds\tau `$ (classical behavior), and $`V^\mu w^\mu `$, when $`ds\tau `$ (fractal behavior). Recalling that $`D=2`$ plays the role of a critical dimension, we stress that, in the asymptotic scaling domain, $`w^\mu (ds/\tau )^{1/2}`$ , in agreement with Ref. , which allows to identify the fractal domain with the quantum one. In the present paper, only these simplest scale laws with a fractal dimension $`D=2`$ are considered.

The above description strictly applies to an individual fractal trajectory. Now, one of the geometric consequences of the non-differentiability and of the subsequent fractal character of space-time itself (not only of the trajectories) is that there is an infinity of fractal geodesics relating any couple of points of this fractal space-time. It has therefore been suggested  that the description of a quantum mechanical particle, including its property of wave-particle duality, could be reduced to the geometric properties of the set of fractal geodesics that corresponds to a given state of this “particle”. In such an interpretation, we do not have to endow the “particle” with internal properties such as mass, spin or charge, since the “particle” is not identified with a point mass which would follow the geodesics, but its internal properties can simply be defined as geometric properties of the fractal geodesics themselves. As a consequence, any measurement is interpreted as a sorting out of the geodesics of which the properties correspond to the resolution scale of the measuring device (as an example, if the “particle” has been observed at a given position with a given resolution, this means that the geodesics which pass through this domain have been selected) .

The transition scale appearing in Eq. (13) yields two distinct behaviors of the system (particle) depending on the resolution at which it is considered. Equation (13) multiplied by $`ds`$ gives the elementary displacement, $`dX^\mu `$, of the system as a sum of two terms

$$dX^\mu =dx^\mu +d\xi ^\mu ,$$
(14)
$`d\xi ^\mu `$ representing the “fractal part” and $`dx^\mu `$ the “classical part”, each term being defined as

$$dx^\mu =v^\mu ds,$$
(15)
$$d\xi ^\mu =a^\mu \sqrt{2𝒟}(ds^2)^{1/2D},$$
(16)
which becomes, for $`D=2`$

$$d\xi ^\mu =a^\mu \sqrt{2𝒟}ds^{1/2},$$
(17)
with $`2𝒟=\tau `$. We note, from Eqs. (14) to (17), that $`dx^\mu `$ scales as $`ds`$, while $`d\xi ^\mu `$ scales as $`ds^{\frac{1}{2}}`$. Therefore, the behavior of the system is dominated by the $`d\xi ^\mu `$ term in the non-differentiable “fractal” domain (below the transition scale), and by the $`dx^\mu `$ one in the differentiable “classical” domain (above the transition scale).

Now, the Klein-Gordon and Dirac equations give results applying to measurements performed on quantum objects, but realised with classical devices, in the differentiable “classical” domain. The microphysical scale at which the physical systems under study are considered induces the sorting out of a bundle of geodesics, corresponding to the scale of the systems, while the measurement process implies a smoothing out of the geodesic bundle coupled to a transition from the non-differentiable “fractal” to the differentiable “classical” domain. We therefore define an operator $`\overline{𝒞\mathrm{}}`$, which we apply to the fractal variables or functions each time we are drawn to the “classical” domain where the $`ds`$ behavior dominates. The effect of $`\overline{𝒞\mathrm{}}`$ is to extract, from the fractal variables or functions to which it is applied the “classical part”, i.e., the part scaling as $`ds`$.

### C Differential proper time symmetry breaking

Another consequence of the non-differentiable nature of space-time is the breaking of the local differential proper time reflection invariance. The derivative with respect to the proper time $`s`$ of a differentiable function $`f`$ can be written twofold

$$\frac{df}{ds}=\underset{ds0}{lim}\frac{f(s+ds)f(s)}{ds}=\underset{ds0}{lim}\frac{f(s)f(sds)}{ds}.$$
(18)
The two definitions are equivalent in the differentiable case. One passes from one to the other by the transformation $`dsds`$ (local differential proper time reflection invariance), which is therefore an implicit discrete symmetry of differentiable physics. In the non-differentiable situation, both definitions fail, since the limits are no longer defined. In the new framework of scale relativity, the physics is related to the behavior of the function during the “zoom” operation on the proper time resolution $`\delta s`$, identified with the differential element $`ds`$. Two functions $`f_+^{}`$ and $`f_{}^{}`$ are therefore defined as explicit functions of $`s`$ and $`ds`$

$$f_+^{}(s,ds)=\frac{f(s+ds,ds)f(s,ds)}{ds},$$
(19)
$$f_{}^{}(s,ds)=\frac{f(s,ds)f(sds,ds)}{ds}.$$
(20)
When applied to the space-time coordinates, these definitions yield, in the non-differentiable domain, two four-velocities of which the components are fractal functions of the resolution, $`V_+^\mu [x^\mu (s),s,ds]`$ and $`V_{}^\mu [x^\mu (s),s,ds]`$. In order to go back to the “classical” domain and to derive the “classical” velocities appearing in Eq. (15), we smooth out each fractal geodesic in the bundles selected by the zooming process with balls of radius larger than $`\tau `$. This amounts to carry out a transition from the non-differentiable to the differentiable domain and leads to define two “classical” velocity fields now resolution independent: $`V_+^\mu [x^\mu (s),s,ds>\tau ]=\overline{𝒞\mathrm{}}V_+^\mu [x^\mu (s),s,ds]=v_+^\mu [x^\mu (s),s]`$ and $`V_{}^\mu [x^\mu (s),s,ds>\tau ]=\overline{𝒞\mathrm{}}V_{}^\mu [x^\mu (s),s,ds]=v_{}^\mu [x^\mu (s),s]`$. The important new fact appearing here is that, after the transition, there is no longer any reason for these two velocities to be equal. While, in standard mechanics, the concept of velocity was one-valued, we must introduce, for the case of a non-differentiable space, two velocities instead of one, even when going back to the “classical” domain. This two-valuedness of the velocity vector finds its origin in the breaking of the discrete proper time reflection invariance symmetry ($`dsds`$). Therefore, if one reverses the sign of the proper time differential element, $`v_+^\mu `$ becomes $`v_{}^\mu `$. But the choice of a reversed proper time element ($`ds`$) must be as qualified as the initial choice ($`ds`$) for the description of the laws of nature, and we have no way, at this level of the description, to define what is the “ past” and what is the “future”. The only solution to this problem is to consider both the forward $`(ds>0)`$ and backward $`(ds<0)`$ processes on the same footing. Then the information needed to describe the system is doubled with respect to the standard differentiable description.

A natural way of accounting for this doubling consists in using complex numbers and the complex product. This is the origin of the complex nature of the wave function of quantum mechanics, since this wave function can be identified with the exponential of the complex action that is naturally introduced in this framework . The choice of complex numbers to represent the two-valuedness of the velocity is not an arbitrary choice. The use of these numbers is the most natural way to generalize to two dimensions the set of real numbers, and obtain, in the relativistic motion case, the complex Klein-Gordon equation . In the same way, as we shall see in Secs. III and IV, the use of bi-quaternionic numbers (i.e., complex Hamilton’s quaternions) is the most natural way to generalize them to eight dimensions  and thus obtain the bi-quaternionic Klein-Gordon equation, which yields the Dirac equation.

## III Bi-quaternionic Klein-Gordon equation

It has long been known that the Dirac equation naturally proceeds from the Klein-Gordon equation when writen in a quaternionic form , i.e., using the quaternionic formalism, as introduced by Hamilton , and further developed by Conway . However, this remains essentially formal in the standard framework, since there is no first principle reason for which the probability amplitude should be quaternionic. In the present approach, we shall see that its bi-quaternionic nature can be established as a consequence of the nonderivable nature of space-time. Indeed, we propose in the current section to introduce naturally a bi-quaternionic covariant derivative operator, leading to the definition of a bi-quaternionic velocity and wave-function, which we use to derive the Klein-Gordon equation in a bi-quaternionic form. This allows us to obtain, in Sec. IV, the Dirac equation as a mere consequence, since Dirac spinors and bi-quaternions are actually two equivalent representations of an electron.

### A Further symmetry breakings and bi-quaternionic covariant derivative operator

Most of the approach described in Sec. II remains applicable. However, the main new features obtained in the case we now study proceed from a deeper description of the scale formalism, considering the more general case when the peculiar choice of an identification of the differentials and the resolution variables (remember the choice explicited in Sec. II B, where we have set $`ds=\delta s`$) is given up, implying the subsequent breaking of the symmetries:

$`dsdsdx^\mu dx^\mu x^\mu x^\mu `$

In the scaling domain, the four space-time coordinates $`X^\mu (s,ϵ_\mu ,ϵ_s)`$ are now four fractal functions of the proper time $`s`$ and of the resolutions $`ϵ_\mu `$ for the coordinates and $`ϵ_s`$ for the proper time. We consider the case when, for an elementary displacement $`dX^\mu `$ corresponding to a shift $`ds`$ in the curvilinear parameter, the resolutions verify $`ϵ_\mu <dX^\mu `$ and $`ϵ_s<ds`$, which implies that, at a given $`X^\mu `$, a forward shift $`ds`$ of $`s`$ yields a displacement $`dX^\mu `$ of $`X^\mu `$ and a backward shift $`ds`$ produces a displacement $`dX^\mu `$, of which the amplitudes are not necessarily equal (see Fig. 1).

We can therefore apply, to these two different elementary displacements, the canonical decomposition

$`dX^\mu =dx^\mu +d\xi ^\mu ,\mathrm{with}\overline{𝒞\mathrm{}}dX^\mu =dx^\mu =v_{\genfrac{}{}{0pt}{}{+}{\mu }}^\mu ds,`$ (21)
$`d\xi ^\mu =a_+^\mu \sqrt{2𝒟}(ds^2)^{\frac{1}{2D}},\overline{𝒞\mathrm{}}a_+^\mu =0,\overline{𝒞\mathrm{}}(a_+^\mu )^2=1,`$ (22)
and

$`dX^\mu =\delta x^\mu +\delta \xi ^\mu ,\mathrm{with}\overline{𝒞\mathrm{}}dX^\mu =\delta x^\mu =v_\genfrac{}{}{0pt}{}{}{\mu }^\mu ds,`$ (23)
$`\delta \xi ^\mu =a_{}^\mu \sqrt{2𝒟}(ds^2)^{\frac{1}{2D}},\overline{𝒞\mathrm{}}a_{}^\mu =0,\overline{𝒞\mathrm{}}(a_{}^\mu )^2=1.`$ (24)
In the differentiable case, $`dX^\mu =(dX^\mu )`$, and therefore $`v_{\genfrac{}{}{0pt}{}{+}{\mu }}^\mu =v_\genfrac{}{}{0pt}{}{}{\mu }^\mu `$. This is no longer the case in the non-differentiable case, where the local symmetry $`dx^\mu dx^\mu `$ is broken.

Furthermore, we must also consider the breaking of the symmetry $`dsds`$ proceeding from the twofold definition of the derivative with respect to the curvilinear parameter $`s`$. Applied to $`X^\mu `$, considering an elementary displacement $`dX^\mu `$, the classical part extraction process gives two classical forward and backward derivatives $`d/ds_+`$ and $`d/ds_{}`$, which yield in turn two classical velocities, which we denote $`v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{+}{\mu }}^\mu `$. Considering now the same extraction process, applied to an elementary displacement $`dX^\mu `$, the classical forward and backward derivatives $`d/ds_+`$ and $`d/ds_{}`$ allow again to define two classical velocities, denoted $`v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{}{\mu }}^\mu `$. We summarize this result as

$$v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{+}{\mu }}^\mu =\frac{dx^\mu }{ds_\pm },v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{}{\mu }}^\mu =\frac{\delta x^\mu }{ds_\pm }.$$
(25)
Contrary to what happens in the differentiable case, the total derivative with respect to the proper time of a fractal function $`f(x(s),s)`$ with integer fractal dimension contains finite terms up to higher order

$$\frac{df}{ds}=\frac{f}{s}+\frac{f}{x^\mu }\frac{dX^\mu }{ds}+\frac{1}{2}\frac{^2f}{x^\mu x^\nu }\frac{dX^\mu dX^\nu }{ds}+\frac{1}{6}\frac{^3f}{x^\lambda x^\mu x^\nu }\frac{dX^\lambda dX^\mu dX^\nu }{ds}+\mathrm{}$$
(26)
In our case, a finite contribution only proceeds from terms of $`D`$-order, while lesser-order terms yield an infinite contribution and higher-order ones are negligible. Therefore, in the special case of a fractal dimension $`D=2`$, the total derivative writes

$$\frac{df}{ds}=\frac{f}{s}+_\mu f.\frac{dX^\mu }{ds}+\frac{1}{2}\frac{^2f}{x^\mu x^\nu }\frac{dX^\mu dX^\nu }{ds}.$$
(27)
Usually the term $`dX^\mu dX^\nu /ds`$ is infinitesimal, but here its classical part reduces to $`\overline{𝒞\mathrm{}}d\xi ^\mu d\xi ^\nu /ds`$.

We can, at this stage, define several total derivatives with respect to $`s`$ of a fractal function $`f`$. We write them - using a compact straightforward notation with summation over repeated indices - after substituting, in the four-dimensional analog of Eq. (27), the expressions for the derivatives of the $`X^\mu `$ obtained when using Eqs. (22) and (24) with the different expressions of Eq. (25) for the classical velocities

$$\frac{df}{ds}{}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{x}\genfrac{}{}{0pt}{}{\pm }{y}\genfrac{}{}{0pt}{}{\pm }{z}\genfrac{}{}{0pt}{}{\pm }{t}}{}^{}=\frac{f}{s}+(v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu +w_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu )\frac{f}{X^\mu }+a_\pm ^\mu a_\pm ^\nu 𝒟\frac{^2f}{dX^\mu dX^\nu },$$
(28)
$$\mathrm{with}w^\mu =a^\mu \sqrt{2𝒟}(ds^2)^{\frac{1}{2D}\frac{1}{2}}.$$
(29)
Now, when we apply the classical part operator to Eq. (28), using Eq. (29) and the properties of the $`a_\pm ^\mu `$ ’s as stated in Eqs. (22) and (24), the $`w^\mu `$ ’s disappear at the first order, but, at the second order, for the fractal dimension $`D=2`$, the fractal behavior of Eq. (16) writes

$$\overline{𝒞\mathrm{}}w_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu w_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\nu }}^\nu =2𝒟\eta ^{\mu \nu }ds,$$
(30)
the $``$ sign in the right-hand side being the inverse of the s-sign in the left-hand side.

Equation (30) follows from the fact that the $`d\xi ^\mu `$’s are of null classical part and mutually independent. The Minkowski metric component, $`\eta ^{\mu \nu }`$ implies indeed that the classical part of every crossed product $`w_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu w_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\nu }}^\nu `$, with $`\mu \nu `$, is null. This is due to the fact that, even if each term in the product scales as $`ds^{1/2}`$, each behaves as an independent fractal fluctuation around its own classical part. Therefore, when we proceed to the smoothing out of the geodesic bundle during the transition from the fractal to the classical domain, we apply a process which is mathematically (not physically) equivalent to a stochastic “Wiener” process, and also more general, since we do not need any Gaussian distribution assumption. Thus, we can apply to the classical part of the $`w`$ crossed product the property of the product of two independent stochastic variables: i.e., the classical part of the product is the product of the classical parts, namely here zero.

Thanks to Eq. (30), the last term of the classical part of Eq. (28) amounts to a Dalambertian, and we can write

$$\frac{df}{ds}{}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{x}\genfrac{}{}{0pt}{}{\pm }{y}\genfrac{}{}{0pt}{}{\pm }{z}\genfrac{}{}{0pt}{}{\pm }{t}}{}^{}=(\frac{}{s}+v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu _\mu 𝒟^\mu _\mu )f,$$
(31)
where the $``$ sign in the right-hand side is still the inverse of the s-sign. When we apply these derivatives to the position vector $`X^\mu `$, we obtain, as expected,

$$\frac{dX^\mu }{ds}{}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}{}^{}=v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu .$$
(32)
We consider now the four fractal functions $`X^\mu (s,ϵ_\mu ,ϵ_s)`$. At this description level, there is no reason for $`(X^\mu )(s,ϵ_\mu ,ϵ_s)`$ to be everywhere equal to $`(X^\mu )(s,ϵ_\mu ,ϵ_s)`$, owing to a local breaking of the P (for $`\mu =x,y,z`$) and T (for $`\mu =t`$) symmetries. Furthermore, as we have stressed above for $`X^\mu `$, at a given $`X^\mu `$, a forward shift $`ds`$ of $`s`$ yields a displacement $`d(X^\mu )`$ of $`X^\mu `$ and a backward shift $`ds`$ produces a displacement $`d(X^\mu )`$, with no necessarily equal amplitudes. We can therefore apply, to these two different elementary displacements, a decomposition similar to the one described in Eqs. (22) and (24), i.e.,

$`d(X^\mu )=\stackrel{~}{d}x^\mu +\stackrel{~}{d}\xi ^\mu ,\mathrm{with}\overline{𝒞\mathrm{}}d(X^\mu )=\stackrel{~}{d}x^\mu =\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{+}{\mu }}^\mu ds,`$ (33)
$`\stackrel{~}{d}\xi ^\mu =\stackrel{~}{a}_+^\mu \sqrt{2𝒟}(ds^2)^{\frac{1}{2D}},\overline{𝒞\mathrm{}}\stackrel{~}{a}_+^\mu =0,\overline{𝒞\mathrm{}}(\stackrel{~}{a}_+^\mu )^2=1,`$ (34)
and

$`d(X^\mu )=\stackrel{~}{\delta }x^\mu +\stackrel{~}{\delta }\xi ^\mu ,\mathrm{with}\overline{𝒞\mathrm{}}d(X^\mu )=\stackrel{~}{\delta }x^\mu =\stackrel{~}{v}_\genfrac{}{}{0pt}{}{}{\mu }^\mu ds,`$ (35)
$`\stackrel{~}{\delta }\xi ^\mu =\stackrel{~}{a}_{}^\mu \sqrt{2𝒟}(ds^2)^{\frac{1}{2D}},\overline{𝒞\mathrm{}}\stackrel{~}{a}_{}^\mu =0,\overline{𝒞\mathrm{}}(\stackrel{~}{a}_{}^\mu )^2=1,`$ (36)
Then, we jump to the classical parts and we are once more confronted to the breaking of the $`dsds`$ symmetry. Applied to $`X^\mu `$, considering an elementary displacement $`d(X^\mu )`$, the classical forward and backward derivatives $`d/ds_+`$ and $`d/ds_{}`$ give again two “classical” velocities, which we denote $`\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{+}{\mu }}^\mu `$. After another classical part extraction, considering now an elementary displacement $`d(X^\mu )`$, the classical forward and backward derivatives $`d/ds_+`$ and $`d/ds_{}`$ yield two other “classical” velocities, denoted $`\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{}{\mu }}^\mu `$. We write in short

$$\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{+}{\mu }}^\mu =\frac{\stackrel{~}{d}x^\mu }{ds_\pm },\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{}{\mu }}^\mu =\frac{\stackrel{~}{\delta }x^\mu }{ds_\pm }.$$
(37)
We therefore obtain new different total derivatives with respect to $`s`$ of a fractal function $`f`$, which we write

$$\frac{\stackrel{~}{d}f}{ds}{}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{x}\genfrac{}{}{0pt}{}{\pm }{y}\genfrac{}{}{0pt}{}{\pm }{z}\genfrac{}{}{0pt}{}{\pm }{t}}{}^{}=\frac{f}{s}+(\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu +\stackrel{~}{w}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu )\frac{f}{X^\mu }+\stackrel{~}{a}_\pm ^\mu \stackrel{~}{a}_\pm ^\nu 𝒟\frac{^2f}{dX^\mu dX^\nu },$$
(38)
$$\mathrm{with}\stackrel{~}{w}^\mu =\stackrel{~}{a}^\mu \sqrt{2𝒟}(ds^2)^{\frac{1}{2D}\frac{1}{2}}.$$
(39)
The classical part operator applied to Eq. (38) yields the classical total derivatives

$$\frac{\stackrel{~}{d}f}{ds}{}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{x}\genfrac{}{}{0pt}{}{\pm }{y}\genfrac{}{}{0pt}{}{\pm }{z}\genfrac{}{}{0pt}{}{\pm }{t}}{}^{}=(\frac{}{s}+\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu _\mu 𝒟^\mu _\mu )f.$$
(40)
When we apply these derivatives to the position vector $`X^\mu `$, we obtain, as expected again,

$$\frac{\stackrel{~}{d}X^\mu }{ds}{}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}{}^{}=\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu .$$
(41)
If we consider the simplest peculiar case when the breaking of the symmetry $`dx^\mu dx^\mu `$ is isotropic as regards the four space-time coordinates (i.e., the signs corresponding to the four $`\mu `$ indices are chosen equal), we are left with eight non-degenerate components - four $`v_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu `$ and four $`\stackrel{~}{v}_{\genfrac{}{}{0pt}{}{\pm }{s}\genfrac{}{}{0pt}{}{\pm }{\mu }}^\mu `$ \- which can be used to define a bi-quaternionic four-velocity,

$`𝒱^\mu =`$ $`{\displaystyle \frac{1}{2}}(v_{++}^\mu +\stackrel{~}{v}_{}^\mu ){\displaystyle \frac{i}{2}}(v_{++}^\mu \stackrel{~}{v}_{}^\mu )+\left[{\displaystyle \frac{1}{2}}(v_+^\mu +v_+^\mu ){\displaystyle \frac{i}{2}}(v_+^\mu \stackrel{~}{v}_{++}^\mu )\right]e_1`$ (42)
$`+`$ $`\left[{\displaystyle \frac{1}{2}}(v_{}^\mu +\stackrel{~}{v}_+^\mu ){\displaystyle \frac{i}{2}}(v_{}^\mu \stackrel{~}{v}_+^\mu )\right]e_2+\left[{\displaystyle \frac{1}{2}}(v_+^\mu +\stackrel{~}{v}_{++}^\mu ){\displaystyle \frac{i}{2}}(\stackrel{~}{v}_+^\mu +\stackrel{~}{v}_+^\mu )\right]e_3.`$ (43)
with $`e_i(i=1,2,3)`$ denoting Hamilton’s imaginary units satisfying the associative but non-commutative algebra

$$e_ie_j=\delta _{ij}+\underset{k=1}{\overset{3}{}}ϵ_{ijk}e_k,i,j=1,2,3$$
(44)
where $`ϵ_{ijk}`$ is the usual completely antisymmetric three-index tensor with $`ϵ_{123}=1`$.

The freedom in the choice of the actual expression for $`𝒱^\mu `$ will be discussed later. It is constrained fy the following requirements: at the limit when $`ϵ_\mu dX^\mu `$ and $`ϵ_sds`$, every $`e_i`$-term in Eq. (43) goes to zero, and, as $`\stackrel{~}{v}_{}^\mu =v_+^\mu `$ in this limit, one recovers the complex velocity of the scale-relativistic Schrödinger equation , $`𝒱^\mu =[v_{++}^\mu +v_+^\mu i(v_{++}^\mu v_+^\mu )]/2`$; at the classical limit, every term in this equation vanishes, save the real term, and the velocity becomes classical in its usual meaning, i.e., real: $`𝒱^\mu =v_{++}^\mu `$.

The bi-quaternionic velocity thus defined corresponds to a bi-quaternionic derivative operator $`{}_{}{}^{}d\text{´}/ds`$, similarly defined, and yielding, when applied to the position vector $`X^\mu `$, the corresponding velocity. For instance, the derivative operator attached to the velocity in Eq. (43) writes

$`{\displaystyle \frac{{}_{}{}^{}d\text{´}}{ds}}={\displaystyle \frac{1}{2}}({\displaystyle \frac{d}{ds}}{}_{++}{}^{}+{\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{}{}^{}){\displaystyle \frac{i}{2}}({\displaystyle \frac{d}{ds}}{}_{++}{}^{}{\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{}{}^{})+[{\displaystyle \frac{1}{2}}({\displaystyle \frac{d}{ds}}{}_{+}{}^{}+{\displaystyle \frac{d}{ds}}{}_{+}{}^{}){\displaystyle \frac{i}{2}}({\displaystyle \frac{d}{ds}}{}_{+}{}^{}{\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{++}{}^{})]e_1`$ (45)
$`+[{\displaystyle \frac{1}{2}}({\displaystyle \frac{d}{ds}}{}_{}{}^{}+{\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{+}{}^{}){\displaystyle \frac{i}{2}}({\displaystyle \frac{d}{ds}}{}_{}{}^{}{\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{+}{}^{})]e_2+[{\displaystyle \frac{1}{2}}({\displaystyle \frac{d}{ds}}{}_{+}{}^{}+{\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{++}{}^{}){\displaystyle \frac{i}{2}}({\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{+}{}^{}+{\displaystyle \frac{\stackrel{~}{d}}{ds}}{}_{+}{}^{})]e_3.`$ (46)
(47)
Substituting Eqs. (31) and (40) into Eq. (47), we obtain the expression for the bi-quaternionic proper-time derivative operator

$$\frac{{}_{}{}^{}d\text{´}}{ds}=[1+e_1+e_2+(1i)e_3]\frac{}{s}+𝒱^\mu _\mu +i𝒟^\mu _\mu ,$$
(48)
the $`+`$ sign in front of the Dalambertian proceeding from the choice of the metric signature $`(+,,,)`$. We keep here, for generality, the $`/s`$ term, stressing that it is actually of no use, since the various physical functions are not explicitly depending on $`s`$. It is easy to check that this operator, applied to the position vector $`X^\mu `$, gives back the bi-quaternionic velocity $`𝒱^\mu `$ of Eq. (43).

It is worth noting that the expression we have written for $`𝒱^\mu `$ in Eq. (43) is one among the various choices we could have retained to define the bi-quaternionic velocity. The main constraint limiting this choice is the recovery of the complex and real velocities at the non-relativistic motion and classical limits. We also choose $`𝒱^\mu `$ such as to obtain the third term in the right-hand side of Eq. (48) under the form of a purely imaginary Dalembertian, which allows to recover an integrable equation of motion. To any bi-quaternionic velocity satisfying both prescriptions corresponds a bi-quaternionic derivative operator $`{}_{}{}^{}d\text{´}/ds`$, similarly defined, and yielding this velocity when applied to the position vector $`X^\mu `$. But, whatever the definition retained, the derivative operator keeps the same form in terms of the bi-quaternionic velocity $`𝒱^\mu `$, as given by Eq. (54). Therefore, the different choices allowed for its definition merely correspond to different mathematical representations leading to the same physical result.

### B Bi-quaternionic wave-function

Since $`𝒱^\mu `$ is bi-quaternionic, the lagrange function is also bi-quaternionic and, therefore, the same is true of the action.

A generalized equivalence principle, as well as a strong covariance principle, allows us to write the equation of motion under a free-motion form, i.e., the differential geodesic equation

$$\frac{{}_{}{}^{}d\text{´}𝒱_\mu }{ds}=0,$$
(49)
where $`𝒱_\mu `$ is the bi-quaternionic four-velocity, e.g., the covariant counterpart of $`𝒱^\mu `$ defined in Eq. (43).

The elementary variation of the action, considered as a functional of the coordinates, keeps the usual form

$$\delta 𝒮=mc𝒱_\mu \delta x^\mu .$$
(50)
We thus obtain the bi-quaternionic four-momentum, as

$$𝒫_\mu =mc𝒱_\mu =_\mu 𝒮.$$
(51)
We are now able to introduce the wave function. We define it as a re-expression of the bi-quaternionic action by

$$\psi ^1_\mu \psi =\frac{i}{cS_0}_\mu 𝒮,$$
(52)
using, in the left-hand side, the quaternionic product derived from Eq. (51), as

$$𝒱_\mu =i\frac{S_0}{m}\psi ^1_\mu \psi .$$
(53)
This is the bi-quaternionic generalization of the definition used in the Schrödinger case: $`\psi =e^{iS/S_0}`$. It is worth stressing here that we could choose, for the definition of the wave function in Eq. (52), a commutated expression in the left-hand side, i.e., $`(_\mu \psi )\psi ^1`$ instead of $`\psi ^1_\mu \psi `$. But with this reversed choice, owing to the non-commutativity of the quaternionic product, we could not obtain the motion equation as a vanishing four-gradient, as in Eq. (61). Therefore, we retain the above simplest choice, as yielding an equation which can be integrated. This non-commutativity induced property might be considered as a key for a future understanding of the parity violation, which will not be developed here.

### C Free-particle Klein-Gordon equation

As, in what follows, we only consider $`s`$-stationary functions, i.e., functions which do not explicitly depend on the proper time $`s`$, the derivative operator reduces to

$$\frac{{}_{}{}^{}d\text{´}}{ds}=𝒱^\nu _\nu +i𝒟^\nu _\nu .$$
(54)
Now this expression is independent of the peculiar choice retained for the bi-quaternionic form of the four-velocity, as the representation dependent term in the right-hand side of Eq. (48) has vanished.

The equation of motion, Eq. (49), thus writes

$$\left(𝒱^\nu _\nu +i𝒟^\nu _\nu \right)𝒱_\mu =0.$$
(55)
We replace $`𝒱_\mu `$, respectively $`𝒱^\nu `$, by their expressions given in Eq. (53) and obtain

$$i\frac{S_0}{m}\left(i\frac{S_0}{m}\psi ^1^\nu \psi _\nu +i𝒟^\nu _\nu \right)\left(\psi ^1_\mu \psi \right)=0.$$
(56)
The particular choice $`𝒮_0=2m𝒟`$ <sup>*</sup><sup>*</sup>*At first sight, the meaning of the choice $`𝒮_0=2m𝒟`$ is not a necessary condition from the viewpoint of physics, but merely a simplifying choice as regards the equation form. Indeed, it is only under this particular choice that the fundamental equation of dynamics can be integrated. However, if we do not make this choice, the $`\psi `$ function is a solution of a third order, non linear, complicated equation such that no precise physical meaning can be given to it. We therefore claim that our choice $`𝒮_0=2m𝒟`$ has a profound physical significance, since, in our construction, the meaning of $`\psi `$ is directly related to the fact that it is a solution of the Klein-Gordon and Dirac equations obtained below. allows us to simplify this equation and we get

$$\psi ^1^\nu \psi _\nu (\psi ^1_\mu \psi )+\frac{1}{2}^\nu _\nu (\psi ^1_\mu \psi )=0.$$
(57)
The definition of the inverse of a quaternion

$$\psi \psi ^1=\psi ^1\psi =1,$$
(58)
implies that $`\psi `$ and $`\psi ^1`$ commute. But this is not necessarily the case for $`\psi `$ and $`_\mu \psi ^1`$ nor for $`\psi ^1`$ and $`_\mu \psi `$ and their contravariant counterparts . However, when we derive Eq. (58) with respect to the coordinates, we obtain

$`\psi _\mu \psi ^1`$ $`=`$ $`(_\mu \psi )\psi ^1`$ (59)
$`\psi ^1_\mu \psi `$ $`=`$ $`(_\mu \psi ^1)\psi ,`$ (60)
and identical formulae for the contravariant analogues.

Developing Eq. (57), using Eqs. (60) and the property $`^\nu _\nu _\mu =_\mu ^\nu _\nu `$, we obtain, after some calculations,

$$_\mu [(^\nu _\nu \psi )\psi ^1]=0.$$
(61)
We integrate this four-gradient as

$$(^\nu _\nu \psi )\psi ^1+C=0,$$
(62)
of which we take the right product by $`\psi `$ to obtain

$$^\nu _\nu \psi +C\psi =0.$$
(63)
We therefore recognize the Klein-Gordon equation for a free particle with a mass $`m`$, so that $`m^2c^2/\mathrm{}^2=C`$.

## IV Dirac equation

We now use a long-known property of the quaternionic formalism, which allows to obtain the Dirac equation for a free particle as a mere square root of the Klein-Gordon operator (see, e.g., Refs. ).

We first develop the Klein-Gordon equation as

$$\frac{1}{c^2}\frac{^2\psi }{t^2}=\frac{^2\psi }{x^2}+\frac{^2\psi }{y^2}+\frac{^2\psi }{z^2}\frac{m^2c^2}{\mathrm{}^2}\psi .$$
(64)
Thanks to the property of the quaternionic and complex imaginary units $`e_1^2=e_2^2=e_3^2=i^2=1`$, we can write Eq. (64) under the form

$$\frac{1}{c^2}\frac{^2\psi }{t^2}=e_3^2\frac{^2\psi }{x^2}e_2^2+ie_1^2\frac{^2\psi }{y^2}i+e_3^2\frac{^2\psi }{z^2}e_1^2+i^2\frac{m^2c^2}{\mathrm{}^2}e_3^2\psi e_3^2.$$
(65)
We now take advantage of the anticommutative property of the quaternionic units ($`e_ie_j=e_je_i`$ for $`ij`$) to add to the right-hand side of Eq. (65) six vanishing couples of terms which we rearrange as

$`{\displaystyle \frac{1}{c}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{1}{c}}{\displaystyle \frac{\psi }{t}}\right)`$ $`=`$ $`e_3{\displaystyle \frac{}{x}}\left(e_3{\displaystyle \frac{\psi }{x}}e_2+e_1{\displaystyle \frac{\psi }{y}}i+e_3{\displaystyle \frac{\psi }{z}}e_1i{\displaystyle \frac{mc}{\mathrm{}}}e_3\psi e_3\right)e_2`$ (66)
$`+`$ $`e_1{\displaystyle \frac{}{y}}\left(e_3{\displaystyle \frac{\psi }{x}}e_2+e_1{\displaystyle \frac{\psi }{y}}i+e_3{\displaystyle \frac{\psi }{z}}e_1i{\displaystyle \frac{mc}{\mathrm{}}}e_3\psi e_3\right)i`$ (67)
$`+`$ $`e_3{\displaystyle \frac{}{z}}\left(e_3{\displaystyle \frac{\psi }{x}}e_2+e_1{\displaystyle \frac{\psi }{y}}i+e_3{\displaystyle \frac{\psi }{z}}e_1i{\displaystyle \frac{mc}{\mathrm{}}}e_3\psi e_3\right)e_1`$ (68)
$``$ $`i{\displaystyle \frac{mc}{\mathrm{}}}e_3\left(e_3{\displaystyle \frac{\psi }{x}}e_2+e_1{\displaystyle \frac{\psi }{y}}i+e_3{\displaystyle \frac{\psi }{z}}e_1i{\displaystyle \frac{mc}{\mathrm{}}}e_3\psi e_3\right)e_3.`$ (69)
We see that Eq. (69) is obtained by applying twice to the bi-quaternionic wavefunction $`\psi `$ the operator

$$\frac{1}{c}\frac{}{t}=e_3\frac{}{x}e_2+e_1\frac{}{y}i+e_3\frac{}{z}e_1i\frac{mc}{\mathrm{}}e_3()e_3.$$
(70)
The three first Conway matrices $`e_3()e_2`$, $`e_1()i`$ and $`e_3()e_1`$ , figuring in the right-hand side of Eq. (70), can be written in the compact form $`\alpha ^k`$, with

$$\alpha ^k=\left(\begin{array}{cc}0& \sigma _k\\ \sigma _k& 0\end{array}\right),$$

the $`\sigma _k`$’s being the three Pauli matrices, while the fourth Conway matrix

$$e_3()e_3=\left(\begin{array}{cccc}1& 0& \mathrm{\hspace{0.33em}\hspace{0.33em}0}& 0\\ 0& 1& \mathrm{\hspace{0.33em}\hspace{0.33em}0}& 0\\ 0& 0& 1& 0\\ 0& 0& \mathrm{\hspace{0.33em}\hspace{0.33em}0}& 1\end{array}\right)$$

can be recognized as the Dirac $`\beta `$ matrix. We can therefore write Eq. (70) as the non-covariant Dirac equation for a free fermion

$$\frac{1}{c}\frac{\psi }{t}=\alpha ^k\frac{\psi }{x^k}i\frac{mc}{\mathrm{}}\beta \psi .$$
(71)
The covariant form, in the Dirac representation, can be recovered by applying $`ie_3()e_3`$ to Eq. (71).

The isomorphism which can be established between the quaternionic and spinorial algebrae through the multiplication rules applying to the Pauli spin matrices allows us to identify the wave function $`\psi `$ to a Dirac spinor. Spinors and quaternions are both a representation of the SL(2,C) group. See Ref.  for a detailed discussion of the spinorial properties of bi-quaternions.

## V Conclusion

The last of the three most fundamental motion equations (Schrödinger, Klein-Gordon and Dirac) which are merely postulated in standard quantum mechanics has been established here, in the framework of Galilean scale relativity, as a geodesic equation in a fractal space-time. The change from classical to quantum relativistic motion arises from successive symmetry breakings in the fractal geodesic picture.

First, the complex nature of the wave function is the result of the differential (proper) time symmetry breaking, which is the simplest effect arising from the fractal structure of space (space-time). At this stage, Galilean scale relativity with a complex wave function permits the derivation of both the Schrödinger  and Klein-Gordon  equations.

To go on with the description of the elementary properties encountered in the microphysical world, we have considered here further breakings of fundamental discrete symmetries, due to nonderivability, namely the differential coordinate symmetry ($`dx^\mu dx^\mu `$) breaking and the parity and time reversal symmetry breakings. These new breakings provide a four-complex component wave function (i.e., a eight component wave function), of which the most natural mathematical representation is in term of bi-quaternionic numbers . We therefore obtain the spinorial and the particle anti-particle nature of elementary objects which we can describe as Dirac spinors. Here, spin arises from the isomorphism between the quaternionic and spinorial representations, both of which are different representations of the SL(2,C) group. At this stage, the Klein-Gordon equation is written in a bi-quaternionic form which naturally yields the free Dirac equation.

It is worth stressing that these results only proceed from a restricted use of the scale-relativistic potentialities. We have, in the present work, limited our investigations to the induced effects of scale laws on the equations of motion, in the framework of dilation laws exhibiting a Galilean group structure, i.e., a fractal space-time with a constant fractal dimension $`D=2`$. This is only one of the simplest levels at which the scale-relativistic program can be achieved. We have also made a series of simplifying choices, which we have explained and justified all along the derivation procedure. Some of them have been dictated by physical or experimental considerations, but others only correspond to special cases, provisionally retained, so that other possibilities will have to be explored in the future.


text:
# Phase mixing of a three dimensional magnetohydrodynamic pulse

## 1 Introduction

The phase mixing of incompressible Alfvén waves (Heyvaerts & Priest 1983) has been intensively studied as a possible mechanism of coronal heating in open magnetic structures of the polar regions of the Sun. Also, there are recent observational indications that the same mechanism operates in coronal loops Ofman & Aschwanden (2002). In a compressible plasma, phase mixing of linearly polarized plane Alfvén waves may provide enhanced heating through the nonlinear generation of fast magnetoacoustic waves. This phenomenon has been intensively studied by Malara et al. (1996); Nakariakov et al. (1997, 1998); Botha et al. (2000); Tsiklauri et al. (2001, 2002) both analytically and numerically.

There are three main motivations for further development of this work: The first, as was recently shown by Hood et al. (2002a), is that phase mixing of localized single Alfvén pulses results in a slower, power-law damping as opposed to the standard $`\mathrm{exp}(t^3)`$ one for harmonic Alfvén waves (Heyvaerts & Priest 1983). This suggests that localized Alfvénic perturbations can transport energy higher into the corona than harmonic ones. Also, Hood et al. (2002b) showed that this power-law damping is replaced by even faster power-law decay as the number of pulses is increased. Such perturbations could be generated e.g. by transient events such as solar flares, coronal mass ejections, etc (Roussev et al. 2001). The second follows from the recent demonstration by Tsiklauri & Nakariakov (2002) that in the linear regime and in 3D geometry, generation of compressive perturbations significantly modifies the dynamics of the Alfvénic part of a MHD pulse. In fact, while solving the propagation part of the problem, i.e. considering an ideal plasma limit, it has been shown that an initially localized 3D Alfvén pulse interacting with a 1D transverse inhomogeneity evolves to a MHD pulse with a significant compressible component. This 3D MHD pulse is subject to phase mixing. Therefore, as a natural step forward, in this study we include a finite plasma resistivity in order to investigate quantitatively dissipation of a 3D MHD pulse via phase mixing. The third follows from Parker’s (1991) earlier suggestion that in order for the phase mixing mechanism to work, it requires the presence of a third ignorable coordinate. Although it has been justly observed that having an ignorable coordinate seems unlikely in a filamentary corona, and, also, in the absence of an ignorable coordinate, the “Alfvén-type waves” are coupled to the “fast magnetosonic-type modes”, we believe it was not justified to dismiss the importance of phase mixing without a rigorous, fully three dimensional study.

## 2 The model

In our model we use the MHD equations

$$\rho _tV+\rho (V)V=p(4\pi )^1B\times \mathrm{curl}B,$$
(1)
$$_tB=\mathrm{curl}(V\times B)+\eta \mathrm{\Delta }B,$$
(2)
$$_tp+Vp+\gamma pV=(\gamma 1)\eta (\mathrm{curl}B)^2,$$
(3)
where $`B`$ is the magnetic field, $`V`$ is the plasma velocity, $`\rho `$ is the plasma mass density, and $`p`$ is the plasma thermal pressure. In what follows we use $`5/3`$ for the value of $`\gamma `$. Also, here $`\eta `$ denotes plasma resistivity.

We solve equations (1)-(3) in Cartesian coordinates ($`x,y,z`$). Note that we solve a fully 3D problem retaining variation in the $`y`$-direction, i.e. ($`/y0`$). A uniform magnetic field $`B_0`$ is in the $`z`$-direction. The plasma configuration has a one-dimensional inhomogeneity in the equilibrium density $`\rho _0(x)`$ and temperature $`T_0(x)`$, while the unperturbed thermal pressure, $`p_0`$, is taken to be constant everywhere.

## 3 Numerical Results

Evolution of a purely Alfvénic perturbation allows a full analytical treatment and this is presented in Appendix A. However, if an initial Alfvénic perturbation depends on the $`y`$-coordinate, contrary to the case studied in the Appendix A, the linearized version of the system of equations (1)-(3) no longer allows decoupling of Alfvénic and compressive perturbations. In fact, all three modes, Alfvén, fast and slow magnetosonic, are inter-coupled. Thus, they cannot be separated. Because of the complexity of the problem in this case, we resort to direct numerical simulations of Eqs.(1)-(3) using the numerical code Lare3d (Arber et al. 2001). We employ the following initial conditions:

$$V_y|_{t=0}=\{\begin{array}{cc}AC_A(x)\left[1+\mathrm{cos}(10\pi z)\right]\hfill & \\ \times \mathrm{exp}\left[\alpha _y^2\left(y+\frac{1}{\sqrt{2}\alpha _y}\right)^2\right]\hfill & |z|0.1\text{ ,}\hfill \\ 0\hfill & \text{elsewhere}\hfill \end{array}$$
(4)
The rest of the perturbations at $`t=0`$ are set to zero (including $`B_y`$, in this way we guarantee fulfillment of div$`B=0`$, at $`t=0`$). In Eq.(4) $`A`$ is the amplitude of the wave, $`C_A(x)`$ is the local Alfvén speed (see Appendix A) and $`\alpha _y`$ is a free parameter that controls the strength of the gradient in the $`y`$ direction of the initial perturbation. In effect, $`\alpha _y`$ controls the strength of the coupling between the Alfvénic and compressive (fast and slow magnetosonic) components of the MHD pulse. This particular choice of $`y`$-dependence in Eq.(4) is dictated by a desire to have the maximum of the absolute value of $`_yV_y`$, which is responsible for the coupling to the compressive modes, at the middle of the simulation cube, where the background plasma inhomogeneity is the strongest. Since Lare3d is a fully non-linear numerical code, we simply set small (linear) initial amplitude of the perturbation (given by Eq.(4)) as $`A=0.0001`$, while the output results are normalized to $`A`$. The simulation cube size was set by the limits $`11.0x,y,z11.0`$. Boundary conditions used in all our simulations are zero-gradient in all three spatial dimensions. In fact, our strategy was to forbid any generated wave fronts to reach the boundaries to avoid contamination of the solution. This of course severely restricts simulation end time. However, with the computational resources at our disposal, it was still possible to reach the developed stage of phase mixing. In all our numerical simulations we fix the plasma $`\beta `$ at $`10^3`$ (which is a plausible value for solar corona) and dimensionless resistivity (see Appendix A), $`\eta `$, at $`5\times 10^4`$ (which is the same value used by Hood et al. (2002a)) unless stated otherwise.

In Fig. (1) we track the evolution of the maximal value (over the whole simulation cube at a given time) of $`B_y`$ and $`V_y`$ at the middle of the density inhomogeneity at $`x=0`$, where phase mixing significantly enhances the dissipation. Note that $`max(|V_y(x=0,y,z,t)|)`$ at $`t=0`$ is not unity as it is normalized to $`A`$, not to $`max(|V_y(x=0,y,z,t=0)|)`$. We gather from this graph that although $`B_y`$ was initially zero, it is rapidly generated (on the time scale of $`0.05`$, which is the ratio of length-scale of the pulse $`1/10\pi `$ to $`C_A(0)=0.6`$), and subsequently it decays due to phase mixing. Most of our numerical runs were done with $`200\times 72\times 400`$ resolution (with a non-uniform grid), and we checked convergence of our results by doubling the number of grid points in the $`y`$-direction. Coincidence of the two runs was satisfactory and, in fact, it progressively improves with the increase of the simulation time. Although we were unable to double the resolution in $`x`$ and $`z`$, proof of the convergence of our results comes from the fact (see Fig. (2)) that the analytic result (given by Eq.(16) with an additional pre-factor 0.5 to account for the splitting of initial $`V_y`$ perturbation into two D’Alambert solutions (Tsiklauri & Nakariakov 2002)) and the numerical solution in the case $`\alpha _y=0`$ (when there is a complete decoupling of the Alfvén and the compressive modes) coincide perfectly for large $`t`$s. In fact, at $`t=9.998143`$, the absolute deviation between the two is $`(0.0778360.077785)/0.077785=0.0006`$, which is a remarkably good match for a fully three dimensional simulation. Figure (2) demonstrates that: (A) an increase in the coupling strength parameter, $`\alpha _y`$, results in an overall decrease of the Alfvénic part ($`B_y`$) of the MHD pulse, as expected, and (B) At large times, all lines on the log-log plot tend to become parallel. Therefore, the amplitude of the Alfvénic component of the MHD pulse decays with the same power-law ($`B_yt^{3/2}`$) as in the case of a localized, pure, Alfvén wave (which occurs only when $`\alpha _y=0`$), i.e. coupling to the compressive waves does not affect the decay power-law of the Alfvénic part of the MHD pulse. In order to investigate and quantify point (A) further, in Fig. (3) we plot the efficiency of decay of the Alfvénic part of the MHD pulse, which we define as $`max(B_y(x=0,y,z)/max(B_y(x=0,y,z,\alpha _y=0)`$, as a function of coupling to the compressive modes strength parameter, $`\alpha _y`$ at $`t=10`$. Open diamonds represent numerical simulation results, while a straight line is the best fit, $`1.00.25403\alpha _y`$. We conclude that the efficiency of decay depends linearly upon the degree of localization of the pulse in the homogeneous transverse direction, i.e. the coupling strength parameter, $`\alpha _y`$. As far as point (B) is concerned, apart from the visual manifestation that at large times all lines tend to become parallel to each other, an actual fitting procedure also confirms that, indeed, amplitude of the Alfvénic part of the MHD pulse decays as a power-law ($`t^{3/2}`$) within a few percent error to the power-law index. This is not a trivial result. From this we reach the conclusion that the dynamics (phase mixing) of weakly non-plane Alfvén waves can still be qualitatively understood in terms of the previous 2.5D models. A certain level of understanding of this profound conclusion can be achieved by analyzing a linearized version of the system of equations (1)-(3) for $`\beta =0`$. In this case, coupling between the Alfvén ($`B_y`$) and fast magnetosonic mode ($`B_x`$) (the slow mode is absent in this approximation) is described by the following system of equations:

$`\left[_{tt}^2C_A^2(x)(_{xx}^2+_{zz}^2)\right]B_x=`$
$`\eta _t_{xx}^2B_x+`$ $`C_A^2(x)_{xy}^2B_y`$ (5)
$`\left[_{tt}^2C_A^2(x)(_{yy}^2+_{zz}^2)\right]B_y=`$
$`\eta _t_{xx}^2B_y+`$ $`C_A^2(x)_{xy}^2B_x`$ (6)
From Eqs.(5) and (6) it follows that the terms responsible for the coupling between the Alfvén and fast magnetosonic mode are $`C_A^2B_{x,y}/(L_xL_y)`$ (where $`L_{x,y}`$ are spatial scales in $`x`$ and $`y`$), while the terms responsible for the dissipation are $`\eta B_{x,y}/(t_AL_x^2)`$ (where $`t_A`$ is a typical Alfvén time scale). At large $`t`$s, because of the action of phase mixing in the region of density inhomogeneity, dissipative terms ($`L_x^2`$) are much greater than the coupling terms ($`L_x^1`$). Therefore, for $`t1`$ the dynamics of $`B_y`$ is governed by the dissipation rather than its coupling to the compressive fast mode. This, serves as an explanation why, at large times, all lines in Fig. (2) tend to become parallel, i.e. the amplitude of Alfvénic component of the MHD pulse decays as the same power-law ($`B_yt^{3/2}`$) as in the case of localized, pure, Alfvén wave.

Also, we investigated the dependence of the Alfvénic part of the MHD pulse, which we quantify by $`max(B_y(x=0,y,z,t=10)`$, upon resistivity, $`\eta `$. The results are presented in Fig. 4. On the one hand, an accurate analytic solution (16) for $`B_y`$, which scales as $`1/\sqrt{\eta }`$, was obtained for the case of $`\alpha _y=0`$. On the other hand, point (B), mentioned above, indicates that when $`\alpha _y0`$, the solution (16) should be altered by a further pre-factor that can be determined e.g. from Fig. (3) (in this case, as $`\alpha _y=1`$, it is 0.7188). Therefore, we can see that this way constructed analytic solution gives a tolerably good fit in Fig. (4). The deviation in the small $`\eta `$ end can be easily explained by the fact that the solution (16) was derived in the developed phase mixing approximation, which for such small values of $`\eta `$ would require numerical runs with end-simulation times much larger than $`t=10`$ (the latter is not possible with current computational facilities, as it would require further increase of the simulation box size). Thus, we conclude that the coupling to its compressive part does not alter the $`1/\sqrt{\eta }`$ behavior of the Alfvénic part of a compressible MHD pulse.

In order to fully understand the dynamics of the Alfvénic part of the MHD pulse under our initial conditions as well as to appreciate fully the extent of the efficiency of the phase mixing for the coronal plasma heating, in Fig. (5) we show the evolution of the absolute value of the $`z`$-profile in the middle of the plasma inhomogeneity of $`V_y`$ with time. We gather from this graph that the initial, very localized, profile splits into two D’Alambert solutions propagating in opposite directions along the field (cf. Tsiklauri & Nakariakov (2002)), which then dissipate rapidly towards a diffuse, Gaussian shape.

Finally, we would like to reiterate the importance of plasma compressibility on the MHD wave dynamics in the context of wave-based theories of coronal heating and MHD turbulence (cf. Tsiklauri & Nakariakov (2002)). Our present numerical simulations confirm that when $`\alpha _y0`$, coupling to the Alfvénic part of the MHD pulse quite efficiently generates a density perturbation (which was initially zero) on a sizeable fraction of the initial Alfvén amplitude.

## 4 Conclusions

In summary, we have established that:

* Efficiency of the decay of the Alfvénic part of a compressible MHD pulse is related linearly to the degree of localization of the pulse in the homogeneous transverse direction, i.e. the coupling strength parameter $`\alpha _y`$.
* In the developed stage of phase mixing (for large times), coupling to its compressive part does not alter the $`1/\sqrt{\eta }`$ behavior of the Alfvénic part of a compressible MHD pulse.
* The dynamics of Alfvén waves can still be qualitatively understood in terms of the previous 2.5D models, e.g., Hood et al. (2002a)’s law ($`B_yt^{3/2}`$) for spatially localized, single pulses, survives in the case of 3D compressive MHD waves (in the absence of an ignorable coordinate).
* All of the above implies that phase mixing remains a relevant paradigm for the coronal heating applications in realistic 3D geometry and compressive plasma.

It is useful to give estimates of typical scales of the problem. We choose a typical (normalization) length scale of the transverse inhomogeneity L=10 Mm and Alfvén speed $`C_A=1000`$ km s-1, which are appropriate for example for coronal plumes. The coronal value for shear (Braginskii) viscosity, by which Alfvén waves damp, is estimated to be about $`\eta =1`$ m2 s-1 (e.g. Hood et al. (2002a), which translates into $`10^{10}`$ in the dimensionless units). However, as some recent studies suggest (Nakariakov et al. 1999; Ofman & Aschwanden 2002) the actual value of the shear viscosity is possibly much larger due to micro-turbulence, e.g. $`10^{9.2\pm 3.5}`$ m2 s-1 (Ofman & Aschwanden 2002). Effectively, this implies that the Alfvén waves damp by bulk Braginskii viscosity rather than the shear one. The latter estimate is based on the observations of standing (kink) Alfvén waves in the coronal loops. Thus, its use in our estimates is justified provided the propagating Alfvén waves, studied here, are affected by the same dissipation mechanism as the standing ones. Using Eq.(15) (for $`x=0`$) we can estimate the characteristic time of damping of the amplitude of an Alfvén wave (when it reduces, say, $`e`$-times) due to phase mixing. Thus, we obtain approximately $`2745.22`$ s for $`\eta =1`$ m2 s-1 and $`0.1634.56`$ s for $`\eta =10^{9.2\pm 3.5}`$ m2 s-1. This translates into 1647.13 Mm (the Alfvén speed in the middle of inhomogeneity is 600 km s-1) for the damping length, for $`\eta =1`$ m2 s-1, and $`0.1020.74`$ Mm for $`\eta =10^{9.2\pm 3.5}`$ m2 s-1. A similar analysis for the harmonic waves (now using Eq.(13) and substituting $`k=2\pi /\sigma `$) yields a damping length of 609.47 Mm for $`\eta =1`$ m2 s-1 and $`0.047.67`$ Mm for $`\eta =10^{9.2\pm 3.5}`$ m2 s-1. These estimates would be affected if we choose different (probably higher) Alfvén speeds or if the pulse is wider. Also, note that in our numerical runs, the width of the Gaussian Alfvén pulse in the longitudinal ($`z`$) direction was fixed at $`\sigma 0.03`$, i.e. 0.3 Mm, while in the transverse ($`y`$) direction the width, $`1/\sqrt{2}\alpha _y`$, was varied from $`+\mathrm{}`$ down to 0.35, i.e. from $`+\mathrm{}`$ down to 3.5 Mm.

###### Acknowledgements.

DT acknowledges financial support from PPARC. Numerical calculations of this work were performed using the PPARC funded Compaq MHD Cluster at St Andrews. GR acknowledges financial support from a Leverhulme Emeritus Fellowship.

## Appendix A

Here, we present an alternative, more rigorous, way to derive the power-law decay of $`B_y`$, suggested by (Hood et al. 2002a). If the initial Alfvénic perturbation (perturbing $`V_y`$ and $`B_y`$) does not depend on the $`y`$-coordinate, the linearized version of the system of equations (1)-(3) allows for a complete decoupling of Alfvénic and compressive perturbations. Therefore, in order to investigate phase mixing quantitatively, we solve the following initial value problem with zero-gradient boundary conditions:

$$\left[_{tt}^2C_A(x)^2_{zz}^2\right]B_y=\eta _t_{xx}^2B_y,$$
(7)
$$B_y|_{t=0}=f(z)$$
(8)
Here, $`C_A(x)`$ is a dimensionless Alfvén speed. In our numerical simulations we have used the following profile for $`C_A(x)`$: $`C_A(x)=1`$ for $`x<0.5`$, $`C_A(x)=0.6+0.4\mathrm{cos}[\pi (x+0.5)]`$ for $`x[0.5,0.5]`$, and $`C_A(x)=0.2`$ for $`x>0.5`$ (we employ here the same normalization as in Tsiklauri & Nakariakov (2002) and, $`\eta `$ is normalized to the product of reference length and reference Alfvén speed as in Hood et al. (2002a)). Note, that as in the developed stage of phase mixing (for large times) $`_{xx}^2_{yy}^2,_{zz}^2`$, in the forthcoming analytic treatment, we only retain largest (over $`x`$-coordinate) derivatives in the Laplacian (while for generality we still keep the full Laplacian in our numerical simulations, see below). Next, we introduce the following Lagrangian coordinates and slow (dissipation) time scale: $`\overline{x}=x`$, $`\xi =zC_A(x)t`$, $`\tau =\epsilon t`$, with $`\epsilon 1`$. In these variables, the leading term on the left hand side of Eq.(7) is $`2\epsilon C_A(x)_{\tau \xi }^2B_y`$, while on the right hand side the leading term is $`\eta C_A(x)_\xi [C_A^{}(x)^2(\tau /\epsilon )^2_{\xi \xi }^2B_y]`$. Here, prime denotes a derivative over $`x`$. After integration over $`\xi `$ and introduction of yet another auxiliary variable, $`s=\eta C_A^{}(x)^2\tau ^3/(6\epsilon ^3)=\eta C_A^{}(x)^2t^3/6`$, we obtain the following (diffusion) equation for $`B_y`$:

$$_sB_y=_{\xi \xi }^2B_y,$$
(9)
which can be easily integrated:

$$B_y=\frac{1}{2\sqrt{\pi s}}_{\mathrm{}}^+\mathrm{}\mathrm{exp}\left[\frac{(\xi \xi ^{})^2}{4s}\right]B_y(\xi ^{},t=0)𝑑\xi ^{}.$$
(10)
From Eq.(10) it immediately follows that in the developed stage of phase mixing, i.e. when $`t,s+\mathrm{}`$, the solution reduces to

$$B_y=\frac{1}{2\sqrt{\pi s}}_{\mathrm{}}^+\mathrm{}B_y(\xi ^{},t=0)𝑑\xi ^{}.$$
(11)
Therefore, provided the integral

$$_{\mathrm{}}^+\mathrm{}B_y(\xi ^{},t=0)𝑑\xi ^{},$$
(12)
is finite (which is not always so e.g. for the harmonic initial conditions, see below), $`B_y`$ scales as $`1/\sqrt{s}`$ or equivalently $`B_y\eta ^{1/2}t^{3/2}`$.

If we substitute a harmonic wave with initial conditions $`B_y(\xi ^{},t=0)=\mathrm{exp}(ik\xi ^{})`$ into Eq.(10), simple integration yields

$$B_y=e^{\eta C_A^{}(x)^2t^3k^2/6}e^{ik(zC_A(x)t)},$$
(13)
thus, we recover the well-known Heyvaerts & Priest’s solution. For the spatially localized (in $`z`$) initial conditions

$$B_y|_{t=0}=\{\begin{array}{cc}1+\mathrm{cos}(10\pi z)\hfill & z[0.1,0.1],\hfill \\ 0\hfill & \mathrm{elsewhere}\hfill \end{array}$$
(14)
following Hood et al. (2002a) we expand $`B_y(\xi ^{},t=0)`$ in terms of Hermite polynomials $`B_y(\xi ^{},t=0)=_{n=0}^{\mathrm{}}\alpha _ne^{\xi _{}^{}{}_{}{}^{2}/2\sigma ^2}_n(\xi ^{}/\sigma )`$ retaining only the fundamental mode $`n=0`$, i.e. to a fairly good accuracy initial conditions (Eq.(14)) can be approximated $`B_y(\xi ^{},t=0)=\alpha _0e^{\xi _{}^{}{}_{}{}^{2}/2\sigma ^2}`$ (note that $`_0=1`$). Therefore, substituting the latter expression for $`B_y(\xi ^{},t=0)`$ into the integral (10), and evaluating the integral analytically we obtain the solution of Eq.(7) to the leading (fundamental $`n=0`$) order

$`B_y={\displaystyle \frac{\alpha _0}{\sqrt{1+\eta C_A^{}(x)^2t^3/3\sigma ^2}}}`$
$`\times \mathrm{exp}\left[{\displaystyle \frac{(zC_A(x)t)^2}{2(\sigma ^2+\eta C_A^{}(x)^2t^3/3)}}\right],`$ (15)
which coincides with the one obtained by Hood et al. (2002a). Here, $`\sigma =0.033`$ and $`\alpha _0=1/5\sqrt{2\pi }\sigma `$. For large time, $`t`$, Eq.(15) implies that the amplitude of the Alfvén wave (e.g. amplitude of the $`B_y`$) decays as

$$B_y=\frac{1}{5}\left[2\pi \eta C_A^{}(x)^2/3\right]^{1/2}t^{3/2},$$
(16)
i.e., in the case of a localized in $`z`$-coordinate Alfvén pulse, Heyvaerts & Priest’s exponential decay is replaced by the power law $`B_yt^{3/2}`$ (Hood et al. 2002a).

Although, the both derivations in the end lead to the same answer, in our opinion, the one presented here is more rigorous, concise and it gives the general solution (Eq.(10)) and its asymptotic form (Eq.(11)) at once.