text:
# Where Are Be/black-hole Binaries?

## 1 Introduction

X-ray binaries consist of a neutron star (NS) or a black hole (BH) accreting from strong stellar winds or via Roche lobe overflow of the companion star. Depending on the masses of the optical companions, they are conventionally divided into low- and high-mass X-ray binaries. By measuring the radial velocity curve of the non-degenerate donor, the value of the mass function can be determined, which provides a lower limit on the mass of the accreting object. When it is combined with the information of the spectra of the donor and the orbital light curves, actual measurements of the mass can be obtained. In this way, currently the masses of 18 compact objects in X-ray binaries have been found to exceed the the maximum mass possible for a neutron star. These binaries are thought to be black-hole X-ray binaries (BHXBs) (see McClintock & Remillard (2003) and references therein).

Among the 18 BHXBs, 3 are persistently bright X-ray systems and 15 are X-ray novae in which 6 show recurrent outbursts. All of the X-ray novae are low-mass X-ray binaries (LMXBs). The three persistent sources, i.e. Cyg X$`1`$, LMC X$`1`$ and LMC X$`3`$ are high-mass X-ray binaries (HMXBs) with massive O/B companion stars. In the most recent catalogue of high mass X-ray binaries (HMXBs) edited by Liu, Paradijs & van den Heuvel (2000), only 20 out of 130 HMXBs are O/B supergiant systems including the black hole binaries. Though roughly two thirds of HMXBs are Be/X-ray binaries, and X-ray pulsations have been found in about sixty Be/X-ray binary systems (Ziolkowski, 2002), there are no acknowledged Be/BH binaries. The existence of Be/BH binaries becomes an open problem. The solutions may be related to two questions: Do they exist? if yes, how can they be observed?

The first question is connected with the birth-rate estimation of BH/massive X-ray binaries (BMXBs). Using Monte-Carlo simulation, Raguzova and Lipunov (1999) calculated the number and distribution of binary BHs with Be stars. They obtained the expected number of Be/BH binaries to be of order of unity per $`2030`$ Be/NS systems. The Be/BH binaries were found to be highly eccentric with orbital periods lying in a range of 10 days to several years. The actual number of Be/BH binaries could be even higher since in their work, the lower limit of $`50M_{\mathrm{}}`$ for the progenitor of the BH was adopted, based on the mass of the supergiant Wray 977 (Kaper et al., 1995). But current stellar evolution models suggest that stars with mass higher than about $`2025M_{\mathrm{}}`$ will result in BHs(Fryer, 1999; Fryer & Kakogera, 2001). More recently Podsiadlowski et al. (2003) did the binary population synthesis calculation with the mass of the BH progenitors in the range of $`25M_{\mathrm{}}<M_\mathrm{p}<45M_{\mathrm{}}`$, and presented the calculated distribution of the orbital period and donor mass for BH binaries in their Fig.2. It is shown that BHs have large possibility to be formed in relatively narrow systems with orbital periods of less than tens of days. The formation rate of BMXBs estimated from their result is about $`10^5\mathrm{yr}^1`$. Although the massive donor stars are not constrained to be Be stars in their work, the conclusion of the orbital distribution may also be applied to Be/BH binaries, since the fraction of the Be stars to the massive stars is not small (see Coe (2000) and references therein). Nevertheless, theoretical estimation of this fraction is difficult. It is not yet clear how Be stars eject the circumstellar envelope. While rapid rotation plays a part in disc creation, it is not the sole cause (Porter, 1999). Whether the physical properties of the donor star (metallicity, instability in the photosphere induced by binary rotation, and etc.) have effects on producing the dense disk winds is still unknown (Zamanov et al., 2001).

The observation of supergiant B\[e\]/BH system, e.g. XTE J$`0421+560`$/CI Cam (Robinson, Ivans & Welsh, 2002), may provide a further proof of the existence of Be/BH binaries, since theoretically sgB\[e\]/X-ray binaries would have less probability to exist than Be X-ray binaries due to the much larger mass of sgB\[e\] stars. Thus, the absence of Be/black-hole binaries might be due to some observational selection effects. In this work, we try to find the clue of the second question, that is, why we can hardly observe them. We would start from the interaction between a Be star and its compact companion.

Be/X-ray binaries are characterized by their transient nature. Recently, theoretical work on disc truncation in Be/X-ray binaries (Okazaki & Negueruela, 2001; Okazaki et al., 2002) was invoked to explain the X-ray outbursts in Be/X-ray binaries, and has got considerable success in the application to several systems (Okazaki & Negueruela, 2001), especially 4U 0115+63/V635 Cas (Negueruela & Okazaki, 2001) and A0535+26 (Haigh, Coe & Fabregat, 2003). The key idea is as follows: The neutron star exerts a net negative tidal torque on the viscous decretion disk of the Be star, diminishing the action of the viscous torque outside some critical radius, and thus resulting in the truncation of the disk. The disc matter would then accumulate in the outer rings of the disk till overcoming of the truncation by the effects of global one-armed oscillations or disk warping, etc. The following sudden infall of high density disk matter onto the neutron star causes type II X-ray outbursts. During the episodes between the outbursts, the neutron star could hardly be observed because of low mass accretion rate or the propeller effect. On the other hand, if the tidal truncation is not very efficient, and the disc can extend beyond the Roche lobe of the Be star at periastron, the matter could be accreted onto the neutron star during the periastron passage, resulting in (quasi-)periodic type I bursts.

It is interesting to see what would happen when the neutron star in the above picture is replaced by a black hole. We extend the application of this theory to arbitrary systems to find the influence of mass ratio, orbital period, eccentricity and viscosity coefficient on the efficiency of disc truncation. The model and the necessary parameters are briefly described in section 2. The results are shown in section 3. Our discussion and conclusions are made in section 4.

## 2 Model

Negueruela & Okazaki (2001) suggest a scenario for the behavior of Be/X-ray binaries based on long-term multiwavelength monitoring and two theoretical models: (1) the viscous decretion disk of Be stars proposed by Lee, Saio & Osaki (1991), where the mass flow is outward in contrast with that in accretion disk, and angular momentum of the disk is added by viscous torque $`T_{\mathrm{vis}}`$; (2) the tidal truncation model investigated by Artymowicz & Lubow (1994) for the gravitational interaction of an eccentric binary system with circumstellar or circumbinary disks.

Following Negueruela & Okazaki (2001), we consider the picture of the Be/X-ray binary in which the compact star of mass $`M_\mathrm{x}`$ moves around a Be star of mass $`M_{}`$ and radius $`R_{}`$ in an orbit of eccentricity $`e`$, exerting resonant torques $`T_{\mathrm{res}}`$ on the Be decretion disk. As first approximation, the disc is assumed to be coplanner. The criterion for the disc truncation at a given resonance radius $`r_{\mathrm{trunc}}`$ is

$$T_{\mathrm{vis}}+T_{\mathrm{res}}0,$$
(1)
where $`T_{\mathrm{res}}`$ can be easily shown to be dominated by the torque of the inner Lindblad resonance $`T_{\mathrm{res}}_{ml}(T_{\mathrm{ml}})_{\mathrm{ILR}}`$. For near-Keplerian discs, given the expression of viscosity and resonance torques as Eqs.(6) and (8) in Negueruela & Okazaki (2001), the above truncation criterion can be approximated as

$$Re^1=\alpha (\frac{H}{r})^2\frac{\pi a^2}{(GM_{})^2(1+q_\mathrm{x})^{2/3}n^{4/3}_{ml}\frac{m^3}{m1}|\varphi _{\mathrm{ml}}|^2},$$
(2)
where the radius of the inner Lindblad resonance is $`r=((m1)/l)^{2/3}(1+q_\mathrm{x})^{1/3}a`$, $`\alpha `$ is the Shakura-Sunyaev viscosity parameter, $`a`$ is the semmimajor axis of the binary orbit, and $`q_\mathrm{x}M_\mathrm{x}/M_{}`$. The Reynolds number $`Re`$ is determined by the thermal structure of the disk.

In our calculation, an isothermal disc with $`c_\mathrm{s}/V_\mathrm{k}(R_{})4.1\times 10^2(T_\mathrm{d}/T_{\mathrm{eff}})^{1/2}`$ is assumed as in Okazaki & Negueruela (2001), where $`c_\mathrm{s}`$ is the sound speed of the disk, $`V_\mathrm{k}`$ is the Keplerian velocity, and the disk temperature $`T_\mathrm{d}`$ is about $`1/2`$ of the effective temperature $`T_{\mathrm{eff}}`$ of the Be star. For a Keplerian disk, the scale height $`H`$ of the disk is then

$$\frac{H}{r}=\frac{c_\mathrm{s}}{V_\mathrm{k}(R_{})}(\frac{r}{R_{}})^{1/2}.$$
(3)
We adopt $`R_{}/R_{\mathrm{}}(M_{}/M_{\mathrm{}})^{0.8}`$ to get the corresponding radius of the donor star.

The key function in the criterion (2) is the potential component $`\varphi _{\mathrm{ml}}`$, which can be expressed as

$$\varphi _{\mathrm{ml}}=\frac{GM_\mathrm{x}}{a}\frac{1}{\pi }_{0}^{}{}_{}{}^{\pi }𝑑f\frac{(1e^2)^{1/2}}{1+e\mathrm{cos}(f)}\mathrm{cos}(mfn(m1)M(f))b_{1/2}^m(r/r_2(f)),$$
(4)
where

$$b_{1/2}^m(x)\frac{2}{\pi }_0^\pi \frac{\mathrm{cos}m\phi d\phi }{\sqrt{12x\mathrm{cos}\phi +x^2}}$$
(5)
is the Laplace coefficient with argument $`x=r/r_2`$, and $`r_2=a(1e^2)/(1+e\mathrm{cos}(f))`$ is the distance of the compact star from the donor star. To make it convenient to solve the integration in $`\varphi _{\mathrm{ml}}`$, we have used the following relations among the mean anomaly $`M`$, the true anomaly $`f`$ and the eccentric anomaly $`E`$ in an orbital ellipse,

$$dM/df=(1e^2)^{3/2}/(1+e\mathrm{cos}(f))^2$$
(6)
$$M=Ee\mathrm{sin}E$$
(7)
$$\mathrm{sin}E=\frac{r_2\mathrm{sin}f}{a\sqrt{1e^2}}=\frac{\sqrt{1e^2}}{1+e\mathrm{cos}f}\mathrm{sin}f.$$
(8)
From Eqs. (7) and (8), we get

$$M(f)=\{\begin{array}{cc}\mathrm{arcsin}(z(f))ez(f),\text{if}0f<\mathrm{arccos}(e);\hfill & \\ (\pi \mathrm{arcsin}(z(f)))ez(f),\text{if}\mathrm{arccos}(e)f\pi \hfill & \end{array}$$
(9)
where $`z(f)=(1e^2)^{1/2}\mathrm{sin}f/(1+e\mathrm{cos}f)`$. Inserting $`M(f)`$ to equation (4), the numerical integration can be done directly.

In the calculation of the truncation criterion (2), $`l=n(m1)`$ for the inner Lindblad resonance, the parameter of the summation is only $`m`$ then. Because the high order components contribute little, we sum the inner Lindblad resonance torque from $`m=2`$ to the value at which the component is three orders smaller than that of $`m=2`$. The truncation radius deduced from the inner Lindblad radius is then

$$r_{\mathrm{trunc}}=n^{2/3}(1+q_\mathrm{x})^{1/3}a=n^{2/3}(GM_{}/(2\pi /P_{\mathrm{orb}})^2)^{1/3}.$$
(10)
The timescale to open the gap between $`r_{\mathrm{trunc}}`$ and the inner Lagrangian point $`d_{\mathrm{L1}}`$ with $`\mathrm{\Delta }r=d_{\mathrm{L1}}r_{\mathrm{trunc}}`$, is about $`t_{\mathrm{open}}Re_{\mathrm{crit}}(\mathrm{\Delta }r/r_{\mathrm{trunc}})^2P_{\mathrm{orb}}/2\pi `$ as in Artymowicz & Lubow (1994), where $`Re_{\mathrm{crit}}=\alpha _{\mathrm{crit}}^1(H/r)^2`$ denotes the Reynolds number for which the gap is opened from $`r_{\mathrm{trunc}}`$. Since the outflow in the disc is subsonic (Okazaki & Negueruela, 2001), the timescale $`\tau _{\mathrm{drift}}\mathrm{\Delta }r/v_\mathrm{r}`$ of a particle drifting from $`r_{\mathrm{trunc}}`$ to the Roche lobe will be longer than the truncation timescale. Thus, the efficient truncation is defined as (Okazaki & Negueruela, 2001)

$$\gamma (\tau _{\mathrm{drift}}/P_{\mathrm{orb}})_{\mathrm{min}}\mathrm{\Delta }r/(0.1c_\mathrm{s}P_{\mathrm{orb}})>1$$
(11)
where $`v_{\mathrm{r}}^{}{}_{\mathrm{max}}{}^{}0.1c_\mathrm{s}`$ and $`d_{\mathrm{L1}}=(0.5000.227\mathrm{lg}q_\mathrm{x})a(1e)`$ \- the distance of the inner Lagrangian point point from the center of the donor (Frank, King & Raine, 2002) at periastron have been used instead of the Roche radius, since it is the flat disk rather than the star itself expanding to the Roche lobe.

## 3 Result

To study the influence of mass ratio, orbital period, eccentricity and viscosity coefficient on the efficiency of disc truncation, we have calculated the cases with $`M_\mathrm{x}=1.4,10M_{\mathrm{}}`$, $`M_{}=6,15,20M_{\mathrm{}}`$, $`P_{\mathrm{orb}}=5,30,100,250\mathrm{d}`$, $`\alpha =0.03,0.1,0.3`$ and $`e=0.01,0.05,0.1,0.3,0.5,0.7`$. The masses of the Be star and the orbital periods are chosen according to Fig. 2. of Podsiadlowski et al. (2003) and the observations of HMXBs. It is shown in their figure that most BH X-ray binaries are born in orbits with $`P_{\mathrm{orb}}<30\mathrm{d}`$ and $`M_{}<15M_{\mathrm{}}`$. On the other hand, most observed Be/NS binaries have the donor stars with masses of about $`1020M_{\mathrm{}}`$, and orbital periods lying in the range of $`16\mathrm{d}<P_{\mathrm{orb}}<300\mathrm{d}`$. The viscosity coefficient and the eccentricity are adopted for convenient comparison with the result in Table 2 of Okazaki & Negueruela (2001).

The calculated values of the truncation efficiency criterion $`(\tau _{\mathrm{drift}}/P_{\mathrm{orb}})_{\mathrm{min}}`$ are shown against the eccentricity in Figs. 1, 2 and 3 with $`M_{}=6,15,20M_{\mathrm{}}`$ respectively. The data of different $`P_{\mathrm{orb}}`$ are grouped. The general features of the figures can be summarized as follows:
(1) the smaller viscosity parameter of the disk, the more efficient truncation (it can be understood from equation (1) that, given the resonance torques by the compact star, discs of weak viscosity are truncated at small radius);
(2) the narrower systems, the more efficient truncation, since $`\gamma a/P_{\mathrm{orb}}P_{\mathrm{orb}}^{1/3}`$;
(3) for systems with $`e>0.1`$, smaller $`e`$ leads to more efficient truncation, because of the corresponding larger distance to the inner lagrangian point at periastron;
(4) the larger mass of the compact star have a little less effective truncations on the disk of Be star (since $`d_{\mathrm{L1}}`$ of Be star is smaller for more massive companion stars, whereas $`r_{\mathrm{trunc}}`$ has no relation to $`M_\mathrm{x}`$ as shown in equation (10)); discs of more massive Be stars would be truncated more effectively, especially when $`e<0.4`$.

Features (1) and (3) have been found in Okazaki & Negueruela (2001), and the viscous effect has been studied in detail by 3D SPH simulations in Okazaki et al. (2002). However, feature (2) was neglected. It is shown in our figures that the influence of the orbital period on the disk truncation is the most promising one among the parameters. Although feature (4) shows that black holes would not have more effective truncation on the decretion disk of the Be stars than neutron stars, with the help of feature (2), we know that if Be/BH binaries are mostly born in narrow systems, effective truncation would still happen (except for very eccentric systems), so that the accretion would be mainly from the polar wind of the donor star.

## 4 Discussion and Conclusion

We have extended the application of Be disk truncation model proposed by Negueruela & Okazaki (2001) and Okazaki & Negueruela (2001) to general Be binary systems, to study the dependance of truncation efficiency on some basic binary parameters, such as the stellar mass, viscosity coefficient, orbital period and eccentricity. The results not only confirm the trend that low viscosity and small eccentricity would lead to effective Be disk truncation (Okazaki & Negueruela, 2001), but also show that most effective truncation would occur in relatively narrow systems.

The above results can help us to understand the observations of Be binaries. In the context of Be/NS binaries, the observed transient features of the wide ($`P_{\mathrm{orb}}17263\mathrm{d}`$) eccentric ($`e0.10.9`$) systems can be well understood by the disk truncation effect. Furthermore, the results also presents a possible explanation for the absence of Be/BH binaries if they are formed in relatively narrow orbits without high eccentricity, as suggested by Podsiadlowski et al. (2003). The relatively narrow orbit of $`P_{\mathrm{orb}}<30\mathrm{d}`$ is supported by the result of binary population synthesis simulation work of Podsiadlowski et al. (2003) with the assumption of circular orbit. The knowledge of the eccentricity could be traced back to the formation of the BHs in the massive binaries. The BHs were born from the collapse of massive stars in two different ways (see Fryer (1999) and reference therein ). (i) The massive stars directly collapse into black holes. (ii) After supernova explosion, a significant part of the stellar matter falls back onto the proto-neutron star forming a black hole. If there were no kick in the latter case, the eccentricity induced by the supernova can be estimated by $`e=\mathrm{\Delta }M/M_\mathrm{f}`$ (Bhattacharya and van den Heuvel, 1991), where $`\mathrm{\Delta }M`$ is the total mass loss from the binary in the supernova explosion, and $`M_\mathrm{f}`$ is the final mass of the binary. Assuming the mass of the secondary star change little during the process, one can expect less $`\mathrm{\Delta }M`$ and larger $`M_\mathrm{f}`$ in black hole binaries, resulting in smaller eccentricity compared to neutron star binaries. The spatial distribution and the kinematics of most of the black-hole binaries appear to be consistent with the assumption that small asymmetric kicks are imparted to the black holes (White & van Paradijs, 1996), with at least one exception GRO J1655-40, the first black hole system which has evidence of runaway space velocity of $`112\pm 18\mathrm{km}\mathrm{s}^1`$ in an orbit of $`e=0.34\pm 0.05`$ (Mirabel et al., 2002).

Because of very effective tidal truncation on Be disks in the relatively narrow and low eccentric systems, the flow of the matter towards the BH will be effectively blocked during almost the whole orbital cycle. In such episode, the accretion could be mainly from the polar wind of the donor star. As suggested by Waters & van Kerkwijk (1989), the polar wind of the Be star probably resembles that in OB stars. Such low-density high-velocity wind can hardly form accretion disks around the BHs, and the accretion would follow the classical Bondi-Hoyle-Littleton (BHL) approximation. For a typical main-sequence $`15M_{\mathrm{}}`$ main-sequence star in the system of $`P_{\mathrm{orb}}=10\mathrm{d}`$, if the mass loss rate $`\dot{M}_\mathrm{w}=10^9M_{\mathrm{}}\mathrm{yr}^1`$ and the terminal wind velocity $`v_{\mathrm{}}=2000\mathrm{km}\mathrm{s}^1`$ exist, the accreting compact star would have luminosity of about $`L_\mathrm{x}10^{33}\mathrm{erg}\mathrm{s}^1`$. For the Be/BH binaries, the luminosity would be still less, since the polar wind just takes a small fraction of the mass loss from the donor star. The wind plasma accumulates in the outer rings of the decretion disc till one-armed oscillaiton instability - probably responsible for the V/R variability seen in the Be stars - occurs, which causes the almost total disruption of the Be disk. The resulted large mass-infall onto the black hole would lead to very luminous outbursts. However, no aspect of the truncation model implies that the large mass transfer must occur (Okazaki et al., 2002). If the collapsed disk due to the dynamical instability falls back on to the more massive Be star, such narrow and low-eccentricity transients would be almost out of detection even for the most sensitive instruments, so it should be careful to deal with the Be stars unrelated to the X-ray sources.

For systems in which effective truncation occurs, e.g. XTE J$`1543568`$ and 2S $`1553542`$, the burst activities are very rare (see Okazaki et al. (2002) and reference therein). During the long quiescent stage, Be stars could be identified from the Balmer (and sometimes other) line emission and the associated strong infrared excess. In the HMXB catalogue of Liu, Paradijs & van den Heuvel (2000), there are probably 24 Be/X-ray systems where the nature of the compact star and the orbital period are undetermined. They are signed as Be/X-ray binaries because of their highly variable X-ray characteristics analogous to those of the well studied Be/NS transients, or the optical identification of the donor stars. Most of them have been observed once during the burst states, and then disappeared from X-ray detection because of very low luminosity. Among these undetermined Be/X-ray systems, XTE J$`1739302`$ and AX J$`0052.97158`$ (SMC 32) can be excluded now since the donor star of the former has been identified as an O supergiant (Smith et al., 2003) and 167.8 seconds pulsations have been found in the latter (Yokogawa et al., 2001). To look for Be/BH binaries in the left systems, optical observations in the quiescent stage become important, from which we can confirm the nature of the donor star and get information of the orbital period and the velocity curve of the donor. The latter two are useful in estimating the dynamical mass of the compact star - the common way to determine whether the compact star is a black hole, though the measurements would be very difficult in Be/X-ray binaries, because of not only the relatively low velocity of the donor but also the variability of the circumstellar envelope of Be stars (Raguzova and Lipunov, 1999). The accurate locations of some of the undetermined Be/X-ray binaries could be derived with Chandra, which will enable following up observations in the infrared/optical wavelength. Till now, we have found very limited data of these faint Be/X-ray binaries, and there is at least one system GRS $`1736297`$ went undetected even in the Chandra observations (Wilson et al., 2003).

Even if the X-ray outbursts due to the large mass inflow from the erupted disks are observed, it is still difficult to distinguish between black holes and neutron stars, unless pulsations have been found in the latter. Although there is growing evidence that black hole X-ray binaries display radio emission when they are in the low/hard X-ray state (Reig, Kylafis & Giannios, 2003), this is not a decisive criterion. LS I$`+61\mathrm{°}303`$ is an example, whether the compact star is a neutron star (Zamanov, 1995) or a magnetized black hole (Punsly, 1999) is still in debate.

Be stars can be distinguished from normal O/B stars by the special dense disk wind structures when they are on the main sequence. Once they evolve to fill the Roche lobe, they can hardly be distinguished. Due to the high mass transfer rate from the Roche lobe overflow, black holes in HMXBs would be extremely luminous, even well above the Eddington limit $`L_{\mathrm{Edd}}`$ for a stellar-mass black hole, resulted from either the beaming effect (King et al., 2001) and/or being genuinely super-Eddington (Begelman, 2002), and thus become ones of the ultra luminous X-ray sources (ULXs). This can also be seen in fig. 6. of Podsiadlowski et al. (2003).

Finally, we conclude that by investigating the disk truncation efficiency in Be/compact star binary systems, we propose a possible explanation for the absence of the observation of Be/BH binaries. We show that if most of the Be/BH binaries are born in relatively short low eccentric orbit, the interaction of the binary stars will lead to very effective truncation on the decretion disk around the Be star, so that most of the Be/BH systems would appear as very low luminous X-ray sources.

We thank Atsuo Okazaki for discussions on the calculation of the Be disk truncation model, and Dr. X.-L. Luo and Dr. Zhi Xu for their help in preparing the figures. We acknowledge the support from NSFC grant 1007003 and NKBRSF-G19990754 grant.


format: text

text:
# 1 Introduction

## 1 Introduction

In modern relativistic theory of reactions, the S-matrix is parameterized on the basis of the formalism of inhomogeneous Lorentz group \[1-7\]. This is similar to the nonrelativistic theory with the essential modification - the relativistic spin rotation. Since, in this formalism, the spin state of a particle is set in its rest frame, the concrete description of this state depends on the frame from which the Lorentz transformation to the rest frame is performed. This is the essence of the effect of relativistic spin rotation at the transition from one frame to another \[5-7\]. The relativistic spin rotation is a purely kinematical effect conditioned by the setting of the spin of a particle in its rest frame and by the additional rotation of the spatial axes at the successive Lorentz transformations along noncolinear directions, the latter leading to the nontransitivity of the parallelism in the theory of relativity (see  and references therein).

For two free particles, the total spin results from the addition of the particle spins, defined in the respective rest frames, according to the usual angular momentum addition rules. The operator of the total angular momentum is represented by a sum of the operators of total spin and orbital angular momentum, both commuting with the free Hamiltonian. The price for this quasi nonrelativistic description is the lost relation of the orbital angular momentum operator with the usual field coordinate. Instead, this operator is related with the c.m.s. coordinate introduced by Pryce , Newton and Wigner  and, within the Dirac theory - by Foldy and Wouthuysen  (see also papers ).

In the case of a nonzero vector of relative velocity, the angles of the relativistic rotation of the spins of the two particles are different - they depend on the frame from which the Lorentz transformations to the particle rest frames are performed. This leads to the frame dependence of the fractions of the states with different values of the total spin (in particular, the singlet and triplet fractions for a system of two spin-1/2 particles). Thus, due to the frame dependence of the difference of the spin rotation angles, the square of the total spin of a system of two free particles with a nonzero vector of relative velocity is not a Lorentz invariant.

In Sections 2 and 3 a brief review of the consequences of the relativistic spin rotation for one-particle spin states is given; the connection of the effect of the relativistic spin rotation with the kinematical spatial rotations at successive Lorentz transformations is analyzed. In Section 4, the relation between the relativistic spin rotation and the Thomas precession of the internal angular momentum at the motion of a particle along a curvilinear trajectory is discussed. The precession of the spin polarization vector of a relativistic particle in the magnetic and electric fields is considered in Section 5. In Sections 6 and 7 the effect of the relativistic spin rotation for the states of two free particles with a nonzero vector of relative momentum is considered in detail. The results are summarized in Section 8.

## 2 Relativistic spatial rotation

Let $`A`$ be any four-vector set in a frame $`K`$: $`A=\{A_0,𝐀\},A^2=A_0^2𝐀^2`$. The components of this four-vector in another frame $`K^{}`$, having parallel respective spatial axes and moving with velocity $`𝐯`$ and Lorentz factor $`\gamma =(1𝐯^2/c^2)^{1/2}`$ ($`c`$ is the velocity of light) with respect to the frame $`K`$, are given by the Lorentz transformation:

$$𝐀^{}=𝐀\frac{𝐯}{c}\frac{\gamma }{\gamma +1}(A_0+A_0^{}),A_0^{}=\gamma \left(A_0\frac{𝐯}{c}𝐀\right),$$
(1)
$`A^2=A^2`$. We will denote this transformation as $`A^{}=L(𝐯)A`$.111 In accordance with the principle of relativity, the parallelism of the coordinate axes requires that the velocity vector of the frame $`K`$ with respect to the frame $`K^{}`$ is $`𝐯`$. One can check that the inverse transformation $`L^1(𝐯)=L(𝐯)`$, i.e. $`A=L(𝐯)A^{}`$. The above definition of parallelism is not unique in the case when the direction of the frame velocity coincides with one of the coordinate axes. In such a case, the observers in the frames $`K`$ and $`K^{}`$ have to specify the plane containing their relative velocity vector and one of the remaining two axes. This can be easily done, e.g., with the help of a light signal.  Let us further consider another frame $`K_0`$ associated with a physical object moving with the velocity $`𝐯_\mathrm{𝟎}`$ and Lorentz factor $`\gamma _0`$, having parallel spatial axes with respect to the frame $`K`$. The same object moves with the velocity $`𝐯_{\mathrm{𝟎}}^{}{}_{}{}^{}`$ and Lorentz factor $`\gamma _0^{}`$ with respect to the frame $`K^{}`$ \- we denote $`\stackrel{~}{K}_0`$ the associated frame having parallel spatial axes with respect to the frame $`K^{}`$. We will see that, generally, $`\stackrel{~}{K}_0K_0`$: the parallel axes of the frames $`K^{}`$ and $`K`$, $`K_0`$ and $`K`$ do not imply the parallel axes of the frames $`K_0`$ and $`K^{}`$. The axes of all the three frames could be mutually parallel if only their velocities were colinear.

Considering a particle of four-velocity $`u=p/m`$ (p is the particle four-momentum and m is its mass) at rest in the frame $`\stackrel{~}{K}_0`$ or $`K_0`$: $`u^{}=\{c,0,0,0\}`$, we can transform it to the frame $`K^{}`$: $`u^{}=L(𝐯_0^{})u^{}=\{\gamma _0^{}c,\gamma _0^{}𝐯_0^{}\}`$, and then to the frame $`K`$: $`u=L(𝐯)L(𝐯_0^{})u^{}`$. Comparing the result of the successive transformations with that of the direct one: $`u=L(𝐯_0)u^{}=\{\gamma _0c,\gamma _0𝐯_0\}`$, we get the vector form of the law of the relativistic addition of the velocities $`𝐯_0^{}`$ and $`𝐯`$:

$$𝐯_0=\frac{\gamma \gamma _0^{}}{\gamma _0}\left[𝐯_0^{}\frac{1}{\gamma }+𝐯\frac{\gamma _0^{}+\gamma _0}{\gamma _0^{}(\gamma +1)}\right]=\left[\frac{𝐯_0^{}}{\gamma }+𝐯\left(1+\frac{𝐯𝐯_0^{}}{𝐯^2}\frac{\gamma 1}{\gamma }\right)\right]\left(1+\frac{𝐯𝐯_0^{}}{c^2}\right)^1.$$
(2)
The second equality follows from the transformation of the Lorentz factors (see the second equation in Eq. (1)):

$$\gamma _0=\gamma \gamma _0^{}\left(1+\frac{𝐯𝐯_0^{}}{c^2}\right).$$

At the relativistic addition of velocities in the reverse order: $`u=L(𝐯_0^{})L(𝐯)u^{}`$, we obtain

$$\stackrel{~}{𝐯}_0=\frac{\gamma \gamma _0^{}}{\gamma _0}\left[𝐯\frac{1}{\gamma _0^{}}+𝐯_0^{}\frac{\gamma +\gamma _0}{\gamma (\gamma _0^{}+1)}\right]=\left[\frac{𝐯}{\gamma _0^{}}+𝐯_0^{}\left(1+\frac{𝐯𝐯_0^{}}{𝐯_0^2}\frac{\gamma _0^{}1}{\gamma _0^{}}\right)\right]\left(1+\frac{𝐯𝐯_0^{}}{c^2}\right)^1.$$
(3)
For colinear vectors $`𝐯`$ and $`𝐯_0^{}`$, the results of the addition of velocities in the direct and reverse order coincide. In the case of noncolinear velocities $`𝐯`$ and $`𝐯_0^{}`$, we have $`𝐯_0\stackrel{~}{𝐯}_0`$ (but always $`|𝐯_0|=|\stackrel{~}{𝐯}_0|,\gamma _0=\stackrel{~}{\gamma }_0`$). Thus, the symmetry in the addition of velocities, in general, is absent. The successive Lorentz transformations along noncolinear directions are noncommutative; the pure Lorentz transformations do not form a group.

Let us now compare the results of direct ($`KK_0`$) and successive ($`KK^{}\stackrel{~}{K}_0`$) Lorentz transformations of an arbitrary four-vector $`A`$. The direct Lorentz transformation of this four-vector from the frame $`K`$ to the frame $`K_0`$ gives

$$A^{}=L(𝐯_0)A,$$
(4)
where $`𝐯_0`$ is the velocity of the frame $`K_0`$ in the frame $`K`$, determined by the relativistic addition of the velocities $`𝐯_0^{}`$ and $`𝐯`$. The successive Lorentz transformations from the frame $`K`$ to the frame $`K^{}`$ and then from the frame $`K^{}`$ to the frame $`\stackrel{~}{K}_0`$ lead to the result

$$\stackrel{~}{A}^{}=L(𝐯_0^{})L(𝐯)A.$$
(5)
It is easy to show (see below) that four-vectors $`A^{}`$ and $`\stackrel{~}{A}^{}`$ do not coincide when the velocities $`𝐯_0^{}`$ and $`𝐯`$ are not colinear. Their time components are equal to each other, however the three-vector $`𝐀^{}`$ is turned with respect to the three-vector $`\stackrel{~}{𝐀}^{}`$ by some angle $`\omega `$ around the vector $`[𝐯_0^{}𝐯]`$ . Thus,

$$|𝐀^{}|=|\stackrel{~}{𝐀}^{}|,|A_0^{}|=|\stackrel{~}{A}_0^{}|,𝐀^{}=\widehat{D}_z(\omega )\stackrel{~}{𝐀}^{},$$
(6)
where $`\widehat{D}_z(\omega )`$ is the matrix of an active rotation by the angle $`\omega `$ around the axis $`z`$ parallel to the vector $`[𝐯_0^{}𝐯]`$ (the rotation in Eq. (6) recovers the one introduced in Ref.  after the substitutions $`\widehat{D}_z(\omega )\widehat{D}_z^1(\omega )`$ and $`zz`$). This means that spatial axes of the frames $`\stackrel{~}{K}_0`$ and $`K_0`$ are generally not parallel - the transitivity of parallelism is violated in the theory of relativity .

To find out the rotation angle $`\omega `$, let us follow reference  and assume that the four-vector $`A`$ is a four-velocity of the frame $`K`$, i.e. in this frame $`A=\{c,0,0,0\}`$. At the direct transition from the frame $`K`$ to the frame $`K_0`$, we have from Eq. (4):

$$𝐀^{}=\gamma _2𝐯_0$$
(7)
while, at the successive Lorentz transformations $`KK^{}\stackrel{~}{K}_0`$, we get from Eq. (5):

$$\stackrel{~}{𝐀}^{}=\gamma _0\stackrel{~}{𝐯}_0.$$
(8)
Recall that the velocities $`𝐯_0`$ and $`\stackrel{~}{𝐯}_0`$ represent the results of relativistic addition of the velocities $`𝐯_0^{}`$ and $`𝐯`$ in direct and reverse order given in Eqs. (2) and (3), respectively. According to Eq. (6),

$$𝐯_0=\widehat{D}_z(\omega )\stackrel{~}{𝐯}_0.$$
(9)
Consequently, the angle of the spatial rotation $`\omega `$ is equal to the angle between the velocities $`𝐯_0`$ and $`\stackrel{~}{𝐯}_0`$. Therefore,

$$\mathrm{sin}\omega =\frac{|[𝐯_0\stackrel{~}{𝐯}_0]|}{v_0^2},\mathrm{cos}\omega =\frac{𝐯_0\stackrel{~}{𝐯}_0}{v_0^2},$$

where $`v_0=|𝐯_0|=|\stackrel{~}{𝐯}_0|`$. The positive sign of the angle $`\omega `$ ($`0\omega \pi `$) is fixed by the negative projection of the difference vector $`\stackrel{~}{𝐯}_0𝐯_0`$ on the direction of the vector $`𝐯`$ at nonrelativistic velocities. Using Eqs. (2) and (3), one can show  that

$$\mathrm{sin}\omega =\gamma \gamma _0^{}\frac{vv_0^{}}{c^2}\mathrm{sin}\theta \frac{1+\gamma +\gamma _0^{}+\gamma _0}{(1+\gamma )(1+\gamma _0^{})(1+\gamma _0)},$$
(10)
where

$$\gamma _0=\stackrel{~}{\gamma }_0=\gamma \gamma _0^{}\left(1+\frac{vv_0^{}\mathrm{cos}\theta }{c^2}\right),$$

$`v=|𝐯|`$, $`v_0^{}=|𝐯_0^{}|`$, $`\theta `$ is the angle between the vectors $`𝐯_0^{}`$ and $`𝐯`$ ($`0\theta \pi `$). Taking into account the equality

$$(1+\gamma +\gamma _0^{}+\gamma _0)^2=2(1+\gamma )(1+\gamma _0^{})(1+\gamma _0)(\gamma _0^21)(\gamma _0^21)\mathrm{sin}^2\theta ,$$
(11)
one can express the cosine of the spatial rotation angle in the following forms:

$$\mathrm{cos}\omega =\frac{(1+\gamma +\gamma _0^{}+\gamma _0)^2}{(1+\gamma )(1+\gamma _0^{})(1+\gamma _0)}1=1\frac{(\gamma 1)(\gamma _0^{}1)}{1+\gamma _0}\mathrm{sin}^2\theta .$$
(12)
Clearly, in the case of colinear velocity vectors $`𝐯`$ and $`𝐯_0^{}`$ ($`\mathrm{sin}\theta =0`$), the angle of the relativistic spatial rotation is equal to zero.

## 3 Relativistic spin rotation

The vector of the spin polarization in some frame is constructed by the Lorentz transformation of the four-tensor of the angular momentum $`M_{lm}`$, or the spin four-vector $`S_m=(1/2)e_{mlns}M_{lm}u_s`$, from this frame into the rest frame of a particle ($`e_{mnls}`$ is the antisymmetric unit four-tensor of the fourth range, $`u`$ is the particle four-velocity. Thus, the spin state of a particle depends on the concrete frame from which the transformation to the rest frame of this particle is performed. As a result, the relativistic spin rotation is a particular case of the kinematical spatial rotation considered above and, it is described by the same relations. It should be noted that this fact was not indicated in the book , in which the spatial rotation at successive Lorentz transformations was considered in detail. On the other hand, in papers \[5-7\] the effect of the relativistic spin rotation was discussed without the connection with the general kinematical effect of the relativistic spatial rotation.

Let us assume that $`K_0`$ is the particle rest frame and consider the transformations of the spin four-vector $`S`$ satisfying the condition $`Su=0`$. In the frame $`K`$ the spin four-vector is $`S`$, the particle velocity is $`𝐯_0`$, and, in accordance with the Lorentz transformation from the frame $`K`$ to the rest frame $`K_0`$, the spin polarization vector $`𝜻`$ in the frame $`K`$ (better to say - related to the frame $`K`$) is given by the relation

$$S^{}=L(𝐯_0)S=\{0,𝜻\}.$$
(13)
Note that $`𝜻=\widehat{𝐬}`$ represents a mean value of the spin operator $`\widehat{𝐬}`$ in the particle rest frame. It is related to the normalized spin polarization: $`𝐏=𝜻/s`$ ($`|𝐏|1`$), where $`s`$ is the particle spin.

The spin four-vector $`S^{}`$ in the frame $`K^{}`$ is connected with the spin four-vector $`S`$ in the frame $`K`$ by the Lorentz transformation: $`S^{}=L(𝐯)S`$, where $`𝐯`$ is the velocity of the frame $`K^{}`$ with respect to the frame $`K`$. In accordance with Eq. (5), we find the spin polarization vector in the frame $`K^{}`$ as a result of the transition from the frame $`K^{}`$ to the rest frame $`\stackrel{~}{K}_0`$:

$$\stackrel{~}{S}^{}=L(𝐯_0^{})L(𝐯)S=\{0,\stackrel{~}{𝜻}\},$$
(14)
where $`𝐯_0^{}`$ is the particle velocity in the frame $`K^{}`$. According to Eq. (6),

$$𝜻=\widehat{D}_z(\omega )\stackrel{~}{𝜻},$$
(15)
$`|𝜻|=|\stackrel{~}{𝜻}|`$, $`\zeta _z=\stackrel{~}{\zeta }_z`$ and $`z[𝐯_0^{}𝐯]`$. We see that the angle $`\omega `$ of the relativistic spin rotation around the axis $`z`$ can be calculated directly with the help of Eqs. (10) and (12) for the kinematical spatial angle. The relation (10) was just given in the Stapp paper  on the relativistic theory of polarization phenomena.

Let us consider the dependence of the angle of the relativistic spin rotation $`\omega =\omega (\theta )`$ on the angle $`\theta `$ between the velocities $`𝐯_0^{}`$ and $`𝐯`$. Without loss of generality, both angles $`\theta `$ and $`\omega `$ can be considered in the interval $`0,\pi `$. It follows from Eqs. (10) and (12) that $`\omega (0)=\omega (\pi )=0.`$ The maximal value

$$\omega _{\mathrm{max}}=\omega (\theta _\mathrm{m})=\mathrm{arccos}\left[1\frac{2(\gamma 1)(\gamma _0^{}1)}{(\gamma +1)(\gamma _0^{}+1)}\right]$$

is achieved at

$$\theta _\mathrm{m}=\mathrm{arccos}\left(\sqrt{\frac{(\gamma 1)(\gamma _0^{}1)}{(\gamma +1)(\gamma _0^{}+1)}}\right)\frac{\pi }{2}.$$

Note that $`\omega _{\mathrm{max}}=\frac{\pi }{2}`$, $`<\frac{\pi }{2}`$ and $`>\frac{\pi }{2}`$ when $`(\gamma _0^{}1)(\gamma 1)=8`$, $`<8`$ and $`>8`$, respectively.
At nonrelativistic velocities in the frame $`K^{}`$ ($`v_0^{}/c1,\gamma _0^{}1,\gamma _0\gamma `$) the angle of the spin rotation is very small ($`\omega 1`$) :

$$\omega \frac{\gamma }{\gamma +1}\frac{v_0^{}v}{c^2}\mathrm{sin}\theta .$$
(16)
It follows from Eqs. (10) and (12) that in the ultrarelativistic limit, when $`\gamma _0^{}\mathrm{},\gamma _0/\gamma _0^{}\gamma [1+(v/c)\mathrm{cos}\theta ]`$ :

$$\mathrm{sin}\omega \frac{v}{c}\mathrm{sin}\theta \frac{1+\gamma [1+(v/c)\mathrm{cos}\theta ]}{(1+\gamma )[1+(v/c)\mathrm{cos}\theta ]},\mathrm{cos}\omega 1\frac{\gamma 1}{\gamma [1+(v/c)\mathrm{cos}\theta ]}\mathrm{sin}^2\theta .$$
(17)
At $`v/c1,\gamma 1`$ one has

$$\omega _{\mathrm{max}}\pi 2\sqrt{2/\gamma },\theta _\mathrm{m}\pi \sqrt{2/\gamma }.$$

In this case the sharp dependence on the angle $`\theta `$ in the vicinity of $`\theta =\pi `$ with the width $`1/\sqrt{\gamma }`$ takes place . Outside this vicinity the angle $`\omega `$ is equal approximately to the angle $`\theta `$.

The relations (17) are valid exactly for massless particles (photons, neutrinos). In this situation the spin rotation angle coincides with the aberration angle (the angle between the vectors $`𝐯_0^{}`$ and $`𝐯_0`$); then the spin projection of the particle onto the direction of its momentum (the helicity, the degree of the circular polarization of light) is the relativistic invariant .

Indeed, the angle between the particle velocity in the frame $`K^{}`$ and the particle velocity in the frame $`K`$ is determined, in general, by the relations:

$$\mathrm{sin}\alpha =\frac{|[𝐯_0^{}𝐯_0]|}{v_0^{}v_0}=\frac{v}{c}\mathrm{sin}\theta \frac{\gamma }{\gamma +1}\frac{\gamma _0^{}+\gamma _0}{\sqrt{\gamma _0^21}},$$

$$\mathrm{cos}\alpha =\frac{𝐯_0^{}𝐯_0}{v_0^{}v_0}=\frac{\gamma _0^{}\gamma }{\sqrt{\gamma _0^21}}\left[\frac{v_0^{}}{c}+\frac{v}{c}\mathrm{cos}\theta \frac{v_0^{}}{c}\frac{\gamma 1}{\gamma }\mathrm{sin}^2\theta \right].$$
(18)
One can show that always the angle of the relativistic spin rotation $`\omega `$ at the transformation $`K^{}K`$ is less than the angle of the momentum rotation $`\alpha `$ around the same axis parallel to the vector $`[𝐯_0^{}𝐯]`$. Thus, $`\omega <\alpha `$. However, in the limit

$$v_0^{}/c1,\gamma _0^{}1,\gamma _0=\gamma \gamma _0^{}[1+(v/c)\mathrm{cos}\theta ]1$$

the relations (18) pass into the expressions (17) for massless particles, so then $`\omega =\alpha .`$

## 4  Thomas precession and the effect of the relativistic spin rotation

The Thomas precession of the spin polarization vector is the important effect of the relativistic kinematics; it takes place even in the absence of the direct dynamical action on the particle spin . As a typical example of the Thomas precession, one can consider the precession of the spin of a charged particle with zero magnetic moment at its quasi-classical motion in the external magnetic field .

Due to the orthogonality condition $`S(t)u(t)=0`$ for the spin four-vector $`S(t)`$ and particle four-velocity $`u(t)`$ at each instant time $`t`$, the following equality holds:

$$dS(t)u(t)+S(t)du(t)=0.$$
(19)
In the absence of a direct dynamical action on the spin, the increment $`dS(t)`$ can be only a linear combination of the four-vectors $`S(t)`$ and $`u(t)`$. Since, in the considered case, the former contribution is forbidden due to the conservation of $`S(t)^2`$ ($`S(t)dS(t)=0`$), the increment $`dS(t)`$ is proportional to $`u(t)`$, with the coefficient fixed by the general relation (19). The spin four-vectors at the instant times $`t`$ and $`t+dt`$ are then connected by the relation:

$$dS(t)S(t+dt)S(t)=u(t)\left(S(t)du(t)\right).$$
(20)
The spin polarization vector is determined in the rest frame of a particle at each instant time; in so doing the time $`t`$ is given in the laboratory frame $`K`$. The Lorentz transformations from the laboratory frame $`K`$ to the instantaneous rest frames $`K_0(t)`$ and $`K_0(t+dt)`$ lead to the equalities

$$S^{}(t)=\{0,𝜻(t)\}_{K_0(t)}=L(𝐯(t))S(t)_K;$$

$$S^{}(t+dt)=\{0,𝜻(t+dt)\}_{K_0(t+dt)}=L\left(𝐯(t)+d𝐯(t)\right)S(t+dt)_K.$$
(21)
Now we will consider the connection between the spin polarization vectors $`𝜻(t)`$ and $`𝜻(t+dt)`$. Let us note that, according to Eq. (2) for the relativistic addition of velocities, the velocity of the frame $`K_0(t_1)`$ with respect to the frame $`K_0(t)`$, resulting from the addition of the laboratory velocity $`𝐯(t_1)`$ of the frame $`K_0(t_1)`$ and the velocity $`𝐯(t)`$ of the laboratory frame $`K`$ with respect to the frame $`K_0(t)`$, is equal to

$$\mathrm{\Delta }𝐯^{}=\left[\frac{𝐯(t_1)}{\gamma (t_1)}𝐯(t)\left(1\frac{𝐯(t)𝐯(t_1)}{c^2}\frac{\gamma (t)}{\gamma (t)+1}\right)\right]\left(1\frac{𝐯(t)𝐯(t_1)}{c^2}\right)^1.$$
(22)
At the infinitesimal time increment one has

$$t_1=t+dt,𝐯(t_1)=𝐯(t)+d𝐯.$$

Then it follows from Eq. (22) that

$$d𝐯^{}=\gamma (t)d𝐯+𝐯\frac{𝐯d𝐯}{c^2}\frac{\gamma ^3(t)}{\gamma (t)+1}.$$
(23)
According to the second formula (21), the spin polarization vector of a particle at the instant time $`t+dt`$ is a result of the direct Lorentz transformation of the spin four-vector at the instant time $`t+dt`$ from the laboratory frame $`K`$ to the rest frame of a particle $`K_0(t+dt)`$. On the other hand, taking into account Eq. (20), one can write

$$L(𝐯(t))S(t+dt)_K=\{S(t)du,𝜻(t)\}_{K_0(t)}.$$

Applying the infinitesimal Lorentz transformation $`L(d𝐯^{})`$ to both sides of this equation, it is easy to see that, in the first order in $`dt`$, one arrives at the following equality:

$$L(d𝐯^{})L(𝐯)S(t+dt)_K=\{0,𝜻(t)\}_{\stackrel{~}{K}_0(t+dt)}.$$

Thus, the spin polarization vector at the instant time $`t`$ can be presented as a result of the successive Lorentz transformations firstly from the laboratory frame to the frame $`K_0(t)`$ and then, from the frame $`K_0(t)`$ to the rest frame $`\stackrel{~}{K}_0(t+dt)`$. In accordance with Eqs. (6) and (15), we obtain

$$L(𝐯(t)+d𝐯)=\widehat{D}_z(d\omega )L(d𝐯^{})L(𝐯(t));$$

$$𝜻(t+dt)=\widehat{D}_z(d\omega )𝜻(t).$$
(24)
In Eq. (24) $`\widehat{D}_z(d\omega )`$ is the matrix of the infinitesimal spatial rotation by the angle $`d\omega `$ around the instantaneous axis $`z`$ parallel to the vector $`[d𝐯^{}𝐯]`$.

We deduce from Eq. (24) that the spin polarization vector rotation for the infinitesimal time interval due the Thomas precession coincides with the relativistic spin rotation at the transition from the frame $`K_0(t)`$, in which the particle velocity is $`d𝐯^{}`$, to the laboratory frame $`K`$, in which the particle velocity is $`𝐯(t)+d𝐯`$.

To find the angular velocity of the Thomas precession one should perform the following substitutions in Eq. (10) for the angle of relativistic spatial rotation:

$$𝐯_0^{}d𝐯^{},\gamma _0^{}1,𝐯𝐯(t),\gamma \gamma (t),𝐯_0𝐯(t)+d𝐯,\gamma _0\gamma (t).$$

Taking into account Eq. (23) for $`d𝐯^{}`$, we obtain the following expression for the infinitesimal angle of the rotation of the spin polarization vector at the particle motion along the curvilinear trajectory (in the vector form):

$$d𝝎=\frac{\gamma ^2(t)}{\gamma (t)+1}\frac{1}{c^2}[d𝐯𝐯(t)].$$
(25)
Hence the spin polarization vector satisfies the precession equation:

$$\frac{d𝜻(t)}{dt}=[𝛀_{\mathrm{Th}}(t)𝜻(t)]$$
(26)
with the angular velocity vector of the Thomas precession :

$$𝛀_{\mathrm{Th}}(t)\frac{d𝝎(t)}{dt}=\frac{\gamma ^2(t)}{\gamma (t)+1}\frac{1}{c^2}[𝐯(t)\frac{d𝐯(t)}{dt}].$$
(27)
It can be presented also as

$$𝛀_{\mathrm{Th}}(t)=(\gamma (t)1)𝛀_0(t),$$
(28)
where

$$𝛀_0(t)=\left[𝐥(t)\frac{d𝐥(t)}{dt}\right]$$

is the instantaneous angular velocity of the rotation of the particle momentum, $`𝐥(t)`$ is the unit vector along the momentum direction. According to Eq. (28), the direction of the angular velocity of the Thomas precession is opposite to the direction of the angular velocity of the momentum rotation. It is clear that in the case of rectilinear motion the Thomas precession of spin is absent.

At nonrelativistic velocities, Eq. (27) reduces to the formula

$$𝛀_{\mathrm{Th}}(t)=\frac{1}{2c^2}\left[𝐯(t)\frac{d𝐯(t)}{dt}\right],$$
(29)
which describes, in particular, the Thomas precession of the gyroscope axis.

## 5  Spin precession in the electromagnetic field

Let us consider the precession of the spin polarization vector of a relativistic particle of charge $`q`$, mass $`m`$ and spin $`s`$, moving along a quasiclassical trajectory in the external electromagnetic field. Its magnetic moment

$$\mu =\frac{q\mathrm{}}{2mc}gs,$$
(30)
where $`\mathrm{}`$ is the reduced Planck constant and $`g`$ is the gyromagnetic ratio determining the anomalous part of the magnetic moment $`\mu ^{}=\mu (g2)/g.`$ In particular, in the case of an electron or a negative muon: $`q=|e|,s=1/2,g2+e^2/(\pi \mathrm{}c)`$.

The angular velocity of the spin precession in the electromagnetic field, which is contained in the equation

$$\frac{d𝜻(t)}{dt}=[𝛀(t)𝜻(t)],$$

can be represented as a sum of two terms:

$$𝛀(t)=𝛀_{\mathrm{Th}}(t)+𝛀_{\mathrm{dyn}}(t).$$
(31)
Here $`𝛀_{\mathrm{Th}}(t)`$ is the angular velocity of the Thomas precession, which is described by Eq. (28), and $`𝛀_{\mathrm{dyn}}(t)`$ is the angular velocity of the dynamical precession conditioned by the presence of the nonzero magnetic field $`𝐇^{}(t)`$ in the rest frame of a particle. In accordance with the known result of nonrelativistic quantum mechanics,

$$𝛀_{\mathrm{dyn}}(t)=\frac{1}{\gamma (t)}\frac{qg}{2mc}𝐇^{}(t);$$
(32)
the inverse Lorentz factor in the expression (32) appears due to the time contraction in the rest frame in comparison with the laboratory frame ($`dt=\gamma (t)dt^{}`$). According to the Lorentz transformations of the laboratory magnetic and electric fields $`𝐇(t)`$ and $`𝐄(t)`$ at the point of the particle location, we have

$$𝐇^{}(t)=\gamma (t)\left\{𝐇(t)\frac{\gamma (t)1}{\gamma (t)}𝐥(t)(𝐇(t)𝐥(𝐭))+\left[𝐄(t)\frac{𝐯(t)}{c}\right]\right\},$$
(33)
where $`𝐥(t)=𝐯(t)/v(t)`$. The equation of particle motion

$$\frac{d𝐩(t)}{dt}=q\left(𝐄(t)+\frac{1}{c}[𝐯(t)𝐇(t)]\right)$$
(34)
determines the instantaneous angular velocity of the rotation of the particle momentum:

$$𝛀_0(t)=\frac{q}{mc\gamma (t)}\left\{𝐇(t)𝐥(t)(𝐇(t)𝐥(t))+\frac{\gamma ^2(t)}{\gamma ^2(t)1}\left[𝐄(t)\frac{𝐯(t)}{c}\right]\right\}.$$
(35)
As a result, using Eqs. (28), (31)-(33) and (35), we arrive at the formula of Bargmann, Michel and Telegdi for the total angular velocity of the spin precession of the relativistic particle at the presence of the external electromagnetic field :

$$𝛀(t)=\frac{q}{2mc}\left\{\left(g2+\frac{2}{\gamma (t)}\right)𝐇(t)\frac{\gamma (t)1}{\gamma (t)}(g2)𝐥(t)(𝐥(t)𝐇(t))\right\}$$

$$\frac{q}{2mc}\left(g\frac{2\gamma (t)}{\gamma (t)+1}\right)\left[𝐄(t)\frac{𝐯(t)}{c}\right].$$
(36)
### 5.1  Magnetic field

In the case of a constant transverse homogeneous magnetic field, when $`𝐄=0,\mathrm{𝐇𝐥}=0`$, it follows from Eqs. (35) and (36) that

$$𝛀_0=\frac{q}{mc\gamma }𝐇;𝛀=\frac{q}{mc}\left(\frac{g2}{2}+\frac{1}{\gamma }\right).$$
(37)
Then the angular velocity of the spin precession $`𝛀`$ and the angular velocity of the momentum rotation $`𝛀_0`$ at the circular motion in the magnetic field are connected by the proportionality relation

$$𝛀=\left(\frac{g2}{2}\gamma +1\right)𝛀_0.$$
(38)
In so doing the angular velocity of the spin rotation with respect to the momentum rotation is determined by the anomalous magnetic moment:

$$\mathrm{\Delta }𝛀=\gamma \frac{g2}{2}𝛀_0.$$
(39)
The formula (39) forms the theoretical basis for the measurement of the muon anomalous magnetic moment . Indeed, at the Dirac value $`g=2`$ the projection of the spin polarization vector onto the momentum direction would not change in time. Let the initial spin polarization vector of a muon $`𝜻(t_0)`$ be directed along the momentum $`𝐩(t_0)`$. After $`n`$ total turns of the muon in the magnetic field, the angle between the spin polarization vector $`𝜻(t_0+nT)`$ and the momentum $`𝐩(t_0+nT)=𝐩(t_0)`$, where $`T=2\pi /|𝛀_0|`$ is the period of the circular motion, becomes equal to

$$\mathrm{\Delta }\theta =n\pi \gamma (g2).$$
(40)
### 5.2  Electric field

In the absence of the magnetic field ($`𝐇=0`$) it follows from Eqs. (35) and (36) that

$$𝛀_0(t)=\frac{q}{mc}\frac{\gamma (t)}{\gamma ^2(t)1}\left[𝐄(t)\frac{𝐯(t)}{c}\right],𝛀(t)=\frac{q}{2mc}\left(g\frac{2\gamma (t)}{\gamma (t)+1}\right)\left[𝐄(t)\frac{𝐯(t)}{c}\right].$$
(41)
In so doing,

$$𝛀(t)=\left(\frac{1}{2}(g2)\frac{\gamma ^2(t)1}{\gamma (t)}+\frac{\gamma (t)1}{\gamma (t)}\right)𝛀_0(t).$$
(42)
In case of a plane trajectory of a charged particle in the electric field, the vectors $`𝛀(t)`$ and $`𝛀_0(t)`$ have a constant direction along the normal $`𝐧`$ to the plane of the motion. Then the angle of the precession of the spin polarization vector is described by the formula

$$\theta (t)=\underset{0}{\overset{t}{}}\left[\frac{1}{2}(g2)\frac{\gamma ^2(t^{})1}{\gamma (t^{})}+\frac{\gamma (t^{})1}{\gamma (t^{})}\right]\frac{d\theta _0(t^{})}{dt^{}}𝑑t^{},$$
(43)
where $`\theta _0(t)`$ is the angle between the initial momentum of the particle and its momentum at the instant time $`t`$. In case of a negligible change of the particle kinetic energy at the motion, the connection between the angle of the spin precession and the angle of the momentum rotation takes a simple form:

$$\theta =\left[\frac{1}{2}(g2)\frac{\gamma ^21}{\gamma }+\frac{\gamma 1}{\gamma }\right]\theta _0.$$
(44)
In the nonrelativistic case

$$\theta =\frac{1}{2}(g1)\frac{v^2}{c^2}\theta _0.$$
(45)
The relation (43) may be applied to the plane channeling of positively charged particles (protons, positrons) in bent crystals . The distortion of the trajectory of a charged particle at its motion along the bent channel is conditioned by the average electric field in the plane of the channel which is perpendicular to the particle momentum. In so doing, the angle of the particle deviation $`\theta _0`$ in Eq. (43) coincides with the angle of the bend of the crystal and, the axis of the spin rotation is perpendicular to the plane of the bend of the crystal.

### 5.3  Thomas precession and the spin-orbit interaction

The central electric field has the structure

$$𝐄=\frac{1}{r}\frac{dV(r)}{dr}𝐫,$$
(46)
where $`r=|𝐫|,V(r)`$ is the central potential. According to Eqs. (32) (33) and (45) the angular velocity of the dynamical spin precession in the central electric field can be presented in the form

$$𝛀_{\mathrm{dyn}}=\frac{qg}{2m^2c^2\gamma }\frac{1}{r}\frac{dV(r)}{dr}𝐋,$$
(47)
where $`𝐋=[\mathrm{𝐫𝐩}]`$ is the orbital angular momentum. On the other hand, it follows from Eqs. (28), (35) and (45) that the angular velocity of the Thomas precession is given by the expression

$$𝛀_{\mathrm{Th}}=\frac{q}{m^2c^2}\frac{1}{\gamma +1}\frac{1}{r}\frac{dV(r)}{dr}𝐋$$
(48)
As a result, we have:

$$𝛀_{\mathrm{Th}}=\frac{2}{g}\frac{\gamma }{\gamma +1}𝛀_{\mathrm{dyn}}.$$

For the electron, the gyromagnetic ratio $`g2`$ and, in the nonrelativistic approximation ($`\gamma 1`$), the following relation holds:

$$𝛀_{\mathrm{Th}}=\frac{1}{2}𝛀_{\mathrm{dyn}}.$$
(49)
From the point of view of quantum mechanics, the equality of the spin precession with the angular velocity $`𝛀`$ is obtained through the commutation of the spin operator $`\widehat{𝐬}`$ with the interaction Hamiltonian

$$\widehat{H}_{\mathrm{int}}=\mathrm{}\widehat{𝐬}\widehat{𝛀}=\mathrm{}\widehat{𝐬}\widehat{𝛀}_{\mathrm{dyn}}+\mathrm{}\widehat{𝐬}\widehat{𝛀}_{\mathrm{Th}}.$$

Taking into account Eq. (49), we arrive at the following expression for the spin-orbit interaction of electrons in atoms :

$$\widehat{H}_{\mathrm{int}}=\frac{1}{2}\mathrm{}\widehat{𝐬}\widehat{𝛀}_{\mathrm{dyn}}=\frac{1}{4}\frac{|e|\mathrm{}}{m^2c^2}\frac{1}{r}\frac{dV(r)}{dr}\widehat{𝝈}\widehat{𝐋}.$$
(50)
Here $`\widehat{𝝈}=2\widehat{𝐬}`$ is the Pauli vector operator, $`\widehat{𝐋}`$ is the operator of the orbital angular momentum. Eq. (50) contains the additional factor 1/2 due to the contribution of the Thomas precession which reduces twice the effect of the direct interaction of the magnetic moment of moving electron with the electric field in the atom.

The spin-orbit interaction in atoms leads to the fine splitting of atom levels.

## 6 Effect of the relativistic spin rotation on spin correlations in a two-particle system

The spin correlations in two-particle quantum systems were analyzed in detail as a tool allowing one to measure the space–time characteristics of particle production \[18,20-23\], to study the two–particle interaction and the production dynamics (see \[19,20-22\] and references therein) and to verify the consequences of the quantum–mechanical coherence with the help of Bell–type inequalities .

The spin state of the system of two particles in an arbitrary frame is described by the two-particle density matrix, the elements of which, $`\rho _{m_1m_{1}^{}{}_{}{}^{};m_2m_{2}^{}{}_{}{}^{}}^{(1,2)},`$ are given in the representation of the spin projections of the first and second particle in the corresponding rest frames onto the common coordinate axis $`z`$ (see, e.g., ). However, one should take into account the relativistic spin rotation conditioned by the additional rotation of the spatial axes at the successive Lorentz transformations along noncolinear directions \[5-7\]. As a result, the concrete description of a particle spin state depends on the frame from which the transition to the particle rest frame is performed. Particularly, the total spin composition of the two–particle state with a nonzero vector of relative momentum is generally frame–dependent due to different relativistic rotation angles of the two spins at the transition to the frame moving in the direction which is not colinear with the velocity vectors of both particles.

Usually, it is convenient to consider the spin correlations in the center-of-mass system (c.m.s.) of the particle pair. This is natural at the addition of the two–particle total spin and the relative orbital angular momentum into the conserved total angular momentum. In some cases, however, it may be useful to make transition to the laboratory, e.g., in the case when the particle scatterings are used as their spin analyzers .222In principle, this transition is not necessary since one can transform the four–vectors defining the polarization analyzers first to the pair c.m.s. and then to the respective particle rest frames. Denoting $`M_1`$, $`M_2`$ and $`𝐩_1=𝐤`$, $`𝐩_2=𝐤`$ the masses and c.m.s. momenta of the two particles, their respective c.m.s. velocities are $`𝐯_1=c𝐤/\sqrt{𝐤^2+M_1^2c^2}`$, $`𝐯_2=c𝐤/\sqrt{𝐤^2+M_2^2c^2}`$. We denote the corresponding laboratory velocities as $`𝐯_1^{}`$ and $`𝐯_2^{}`$, and - the laboratory velocity of the particle pair as $`𝐯`$. At the Lorentz transformation from the c.m.s. of the particle pair to the laboratory frame with parallel respective spatial axes, the spins of the first and the second particle (in their respective rest frames) rotate in opposite directions around the axis which is parallel to the vector $`[\mathrm{𝐤𝐯}]`$. The angles of the relativistic spin rotation of the first ($`i=1`$) and second ($`i=2`$) particle are found from Eqs. (10) and (12) with the substitutions

$$𝐯_\mathrm{𝟎}^{}𝐯_i,𝐯_\mathrm{𝟎}𝐯_i^{}$$

and similarly for the corresponding Lorentz factors. Thus,

$$\mathrm{sin}\omega _1=\gamma \gamma _1\frac{vv_1}{c^2}\mathrm{sin}\theta \frac{1+\gamma +\gamma _1+\gamma _1^{}}{(1+\gamma )(1+\gamma _1)(1+\gamma _1^{})},$$

$$\mathrm{cos}\omega _1=1\frac{(\gamma 1)(\gamma _11)}{1+\gamma _1^{}}\mathrm{sin}^2\theta ;$$
(51)
$$\mathrm{sin}\omega _2=\gamma \gamma _2\frac{vv_2}{c^2}\mathrm{sin}\theta \frac{1+\gamma +\gamma _2+\gamma _2^{}}{(1+\gamma )(1+\gamma _2)(1+\gamma _2^{})};$$

$$\mathrm{cos}\omega _2=1\frac{(\gamma 1)(\gamma _21)}{1+\gamma _2^{}}\mathrm{sin}^2\theta .$$
(52)
Here

$$\gamma _1^{}=\gamma \gamma _1\left(1+\frac{vv_1}{c^2}\mathrm{cos}\theta \right),\gamma _2^{}=\gamma \gamma _2\left(1\frac{vv_2}{c^2}\mathrm{cos}\theta \right);$$
(53)
$`\theta =\theta _1`$ is the angle between the vectors $`𝐤`$ ($`𝐯_1`$) and $`𝐯`$; $`v_i=|𝐯_i|`$, $`v_i^{}=|𝐯_i^{}|`$, $`v=|𝐯|`$ and $`\gamma _i`$, $`\gamma _i^{}`$, $`\gamma `$ are the corresponding Lorentz factors. In the case of equal–mass particles the relations $`𝐯_1=𝐯_2,\gamma _1=\gamma _2`$ hold (but $`\gamma _1^{}\gamma _2^{}`$ when $`\theta \pi /2`$).

It should be noted that the positive sign of the angle of the spin rotation corresponds to the direction of the nearest rotation from the vector $`𝐤`$ to the vector $`𝐯`$. It is clear that the angle between the vectors $`𝐤`$ ($`𝐯_2`$) and $`𝐯`$ is equal to $`\theta _2=\theta +\pi `$, and $`\mathrm{sin}\theta _2=\mathrm{sin}\theta `$, $`\mathrm{cos}\theta _2=\mathrm{cos}\theta `$. Without loss of generality one can assume that $`0\theta \pi `$. Since the spins of the first and second particle rotate in opposite directions, the corresponding rotation angles have the opposite signs: $`\omega _1>0,\omega _2<0`$.

In the case of the colinearity of the vectors $`𝐤`$ and $`𝐯`$, when $`\theta =0`$ or $`\theta =\pi `$, both the rotation angles are equal to zero. At nonrelativistic velocities $`v_i`$ in the c.m.s. of the particle pair the angles $`\omega _i`$ of the relativistic spin rotation are very small and scale with $`v_i`$ (see Eq. (16) with the substitutions $`v_0^{}v_1`$ and $`v_0^{}v_2`$).

Taking into account the relativistic spin rotation at the transition from the two–particle c.m.s. to the laboratory, the two–particle spin density matrix is transformed as follows:

$$\widehat{\rho }^{(1,2)}=\widehat{D}^{(1)}(\omega _1)\widehat{D}^{(2)}(\omega _2)\widehat{\rho }^{(1,2)}\widehat{D}^{(1)+}(\omega _1)\widehat{D}^{(2)+}(\omega _2),$$
(54)
where

$$\widehat{D}^{(1)}(\omega _1)=\mathrm{exp}(i\omega _1\widehat{𝐬}_1𝐧),\widehat{D}^{(2)}(\omega _2)=\mathrm{exp}(i\omega _2\widehat{𝐬}_2𝐧)$$
(55)
are the matrices of the active space rotations generated by the vector spin operators $`\widehat{𝐬}_1`$ and $`\widehat{𝐬}_2`$, $`𝐧`$ is the unit vector parallel to the direction of the vector $`[\mathrm{𝐤𝐯}]`$.

In the case of two spin–1/2 particles, the two-particle spin density matrix has the structure \[18,20-22\]:

$$\widehat{\rho }^{(1,2)}=\frac{1}{4}[\widehat{I}^{(1)}\widehat{I}^{(2)}+(\widehat{𝝈}^{(1)}𝐏_1)\widehat{I}^{(2)}+\widehat{I}^{(1)}(\widehat{𝝈}^{(2)}𝐏_2)+\underset{l=1}{\overset{3}{}}\underset{k=1}{\overset{3}{}}T_{lk}\widehat{\sigma }_l^{(1)}\widehat{\sigma }_k^{(2)}].$$
(56)
Here $`\widehat{I}`$ is the two-row unit matrix, $`\widehat{𝝈}^{(i)}=2\widehat{𝐬}_i`$ are the Pauli vector operators, $`𝐏_i=\widehat{𝝈}^{(i)}`$ are the polarization vectors ($`i=1,2`$), $`T_{lk}=\widehat{\sigma }_l^{(1)}\widehat{\sigma }_k^{(2)}`$ are the components of the correlation tensor, $`\{1,2,3\}\{x,y,z\}`$. The left and right indexes of the correlation tensor correspond to the rest frames of the first ($`i=1`$) and second ($`i=2`$) particle, respectively. The corresponding probability to select the particles with the polarizations $`𝜶^{(i)}`$ can be obtained by the substitution of the matrices $`\widehat{\sigma }_l^{(i)}`$ in the expression (56) with the corresponding projections $`\alpha _l^{(i)}`$. Particularly, when analyzing the polarization states with the help of particle decays, the vector analyzing power $`𝜶^{(i)}=\alpha _i𝐧_i`$, where $`\alpha _i`$ is the decay asymmetry corresponding to the decay analyzer unit vector $`𝐧_i`$. As a result \[20-22\], the correlation between the decay analyzers is determined by the product of the decay asymmetries and the trace of the spin correlation tensor

$$T=T_{xx}+T_{yy}+T_{zz}.$$

For example, the angular correlation $`𝐧_1𝐧_2=\mathrm{cos}\theta _{12}`$ between the directions of the three–momenta of the decay protons in the respective rest frames of two $`\mathrm{\Lambda }`$-hyperons decaying into the channel $`\mathrm{\Lambda }p+\pi ^{}`$ with the $`P`$-odd asymmetry $`\alpha =0.642`$ is described by the normalized probability density \[19-22\]

$$W(\mathrm{cos}\theta _{12})=\frac{1}{2}\left(1+\alpha ^2\frac{T}{3}\mathrm{cos}\theta _{12}\right).$$
(57)
Clearly, the structure of both Eq. (56) and the corresponding angular distribution of the spin analyzers (e.g., Eq. (57)) does not depend on the system from which the transitions to the particle rest frames are performed. The system dependence manifests only through the relativistic rotations in the successive Lorentz transformations along noncolinear directions. This circumstance was not understood in paper , where the unnecessary condition of nonrelativistic velocities of $`\mathrm{\Lambda }`$-particles was required.

## 7  Transformations of the spin states of two spin–1/2 particles

The matrices of the space rotations due to the transition from the c.m.s. of two free spin-1/2 particles to the laboratory are the following:

$$\widehat{D}^{(1)}(\omega _1)=\mathrm{cos}\frac{\omega _1}{2}+i\widehat{𝝈}^{(1)}𝐧\mathrm{sin}\frac{\omega _1}{2},\widehat{D}^{(2)}(\omega _2)=\mathrm{cos}\frac{\omega _2}{2}+i\widehat{𝝈}^{(2)}𝐧\mathrm{sin}\frac{\omega _2}{2}.$$
(58)
Selecting the $`z`$–axis parallel to the direction of the vector $`𝐧=[\mathrm{𝐤𝐯}]/|[\mathrm{𝐤𝐯}]|`$, and the axes $`x`$ and $`y`$ in the plane perpendicular to this vector, the polarization vectors and the spin correlation tensor transform at the transition to the laboratory in accordance with the active rotations around the z–axis by the angles $`\omega _1`$ and $`\omega _2`$ for the first and second particle, respectively (the components of the polarization vectors transform according to Eq. (15) with the substitutions: $`\stackrel{~}{𝜻}𝐏_i/2`$ and $`𝜻𝐏_i^{}/2`$):

$$P_{i;x}^{}=P_{i;x}\mathrm{cos}\omega _iP_{i;y}\mathrm{sin}\omega _i;P_{i;y}^{}=P_{i;y}\mathrm{cos}\omega _i+P_{i;x}\mathrm{sin}\omega _i;P_{i;z}^{}=P_{i;z}$$

$$T_{xx}^{}=(T_{xx}\mathrm{cos}\omega _1T_{yx}\mathrm{sin}\omega _1)\mathrm{cos}\omega _2(T_{xy}\mathrm{cos}\omega _1T_{yy}\mathrm{sin}\omega _1)\mathrm{sin}\omega _2;$$

$$T_{yy}^{}=(T_{yy}\mathrm{cos}\omega _1+T_{xy}\mathrm{sin}\omega _1)\mathrm{cos}\omega _2+(T_{yx}\mathrm{cos}\omega _1+T_{xx}\mathrm{sin}\omega _1)\mathrm{sin}\omega _2;T_{zz}^{}=T_{zz};$$

$$T_{xy}^{}=(T_{xy}\mathrm{cos}\omega _1T_{yy}\mathrm{sin}\omega _1)\mathrm{cos}\omega _2+(T_{xx}\mathrm{cos}\omega _1T_{yx}\mathrm{sin}\omega _1)\mathrm{sin}\omega _2;$$

$$T_{yx}^{}=(T_{yx}\mathrm{cos}\omega _1+T_{xx}\mathrm{sin}\omega _1)\mathrm{cos}\omega _2(T_{yy}\mathrm{cos}\omega _1+T_{xy}\mathrm{sin}\omega _1)\mathrm{sin}\omega _2;$$

$$T_{xz}^{}=T_{xz}\mathrm{cos}\omega _1T_{yz}\mathrm{sin}\omega _1;T_{zx}^{}=T_{zx}\mathrm{cos}\omega _2T_{zy}\mathrm{sin}\omega _2;$$

$$T_{yz}^{}=T_{yz}\mathrm{cos}\omega _1+T_{xz}\mathrm{sin}\omega _1;T_{zy}^{}=T_{zy}\mathrm{cos}\omega _2+T_{zx}\mathrm{sin}\omega _2.$$
(59)
Particularly, the trace of the spin correlation tensor transforms at the transition to the laboratory as:

$$T^{}=(T_{xx}+T_{yy})\mathrm{cos}(\omega _1\omega _2)+(T_{xy}T_{yx})\mathrm{sin}(\omega _1\omega _2)+T_{zz}$$
(60)
or, in the case of a symmetric tensor, as:

$$T^{}=T2(T_{xx}+T_{yy})\mathrm{sin}^2\frac{\omega _1\omega _2}{2}.$$
(61)
So, the c.m.s. trace $`T`$ in Eq. (57) is substituted by the laboratory one $`T^{}`$ calculated using Eq. (60) or (61) together with Eqs. (10) and (12) for the spin rotation angles.

It was shown \[20-22\] (see also ) that the trace of the correlation tensor of a system of two spin-1/2 particles is the following linear combination of the relative fractions of singlet (the total spin $`S=0`$) and triplet ($`S=1`$) states:

$$T=\widehat{𝝈}^{(1)}\widehat{𝝈}^{(2)}=\rho _t3\rho _s,\rho _t+\rho _s=1.$$
(62)
When we have the pure singlet state of the particle pair in its c.m.s. ($`\rho _s=1,\rho _t=0,T_{lk}=\delta _{lk},T=3`$), the transformation to the laboratory gives

$$T^{}=3+4\mathrm{sin}^2\frac{\omega _1\omega _2}{2}.$$
(63)
It follows from Eqs. (62) and (63) that at the transition to the laboratory the relative fraction of the singlet state decreases in favor of a triplet state:

$$\rho _s^{}=\mathrm{cos}^2\frac{\omega _1\omega _2}{2},\rho _t^{}=\mathrm{sin}^2\frac{\omega _1\omega _2}{2}.$$
(64)
Thus the square of the total spin of two free particles with a nonzero vector of relative velocity is not a relativistic invariant (see ). Introducing the two–particle singlet state:

$$|\psi _{00}=\frac{1}{\sqrt{2}}\left(|+1/2_z^{(1)}|1/2_z^{(2)}|1/2_z^{(1)}|+1/2_z^{(2)}\right)$$
(65)
and the triplet state with the zero projection onto the rotation axis $`z`$:

$$|\psi _{10}=\frac{1}{\sqrt{2}}\left(|+1/2_z^{(1)}|1/2_z^{(2)}+|1/2_z^{(1)}|+1/2_z^{(2)}\right),$$
(66)
the result in Eq. (64) also follows directly from the matrices of space rotations in Eq. (58); the singlet state in the two–particle c.m.s. is transformed into the following superposition of the singlet and triplet states in the laboratory:

$$|\psi _s^{}=\mathrm{cos}\frac{\omega _1\omega _2}{2}|\psi _{00}+i\mathrm{sin}\frac{\omega _1\omega _2}{2}|\psi _{10}.$$
(67)
Similarly, the transformation of the pure triplet state $`|\psi _{10}`$ in the two–particle c.m.s. ($`\rho _s=0,\rho _t=1,T_{zz}=1,T_{xx}=T_{yy}=1,T=1`$) to the laboratory gives

$$T^{}=14\mathrm{sin}^2\frac{\omega _1\omega _2}{2},$$
(68)
the corresponding fractions being

$$\rho _s^{}=\mathrm{sin}^2\frac{\omega _1\omega _2}{2},\rho _t^{}=\mathrm{cos}^2\frac{\omega _1\omega _2}{2},$$
(69)
in accordance with the transformation:

$$\psi _t^{}=\mathrm{cos}\frac{\omega _1\omega _2}{2}|\psi _{10}+i\mathrm{sin}\frac{\omega _1\omega _2}{2}|\psi _{00}.$$
(70)
In the case of the unpolarized triplet in the two–particle c.m.s. ($`\rho _s=0,\rho _t=1,T_{lk}=(1/3)\delta _{lk},T=1`$ ) we have

$$T^{}=1\frac{4}{3}\mathrm{sin}^2\frac{\omega _1\omega _2}{2},$$
(71)
$$\rho _s^{}=\frac{1}{3}\mathrm{sin}^2\frac{\omega _1\omega _2}{2},\rho _t^{}=1\frac{1}{3}\mathrm{sin}^2\frac{\omega _1\omega _2}{2}.$$
(72)
Using Eqs. (10) and (12), it is easy to show that in the case of two spin-1/2 particles with the same masses ($`\gamma _2=\gamma _1,v_2=v_1`$) the measure of the spin mixing,

$$\kappa =\mathrm{sin}^2\frac{\omega _1\omega _2}{2},$$

can be written in the form:

$$\kappa =\frac{(v_1v)^2}{c^2}\mathrm{sin}^2\theta \left[\left(\frac{1}{\gamma }+\frac{1}{\gamma _1}\right)^2+\frac{(v_1v)^2}{c^2}\mathrm{sin}^2\theta \right]^1.$$
(73)
The maximum of the mixing factor $`\kappa `$ corresponds to the angle $`\theta =\pi /2`$. In the ultrarelativistic limit, when $`\gamma _11,\gamma 1`$ and $`\mathrm{sin}\theta \mathrm{max}(1/\gamma ,\mathrm{\hspace{0.17em}1}/\gamma _1)`$, the factor $`\kappa `$ approaches unity. Then the singlet state in the two–particle c.m.s. becomes in the laboratory the triplet state with the zero projection of the total spin onto the spin rotation axis $`z`$ and, vice versa.

## 8 Summary

1. The effect of the relativistic spin rotation at the transition from some frame to another frame is analyzed. The essence of this effect lies in the fact that in the framework of the formalism of inhomogeneous Lorentz group the particle spin state is set in the particle rest frame and, its concrete description depends on the frame from which the Lorentz transformation to the rest frame is performed.

2. The connection between the relativistic spin rotation and the relativistic spatial rotation at the successive Lorentz transformations along noncolinear directions is considered. It is shown that the relativistic spin rotation is a purely kinematical effect; the rotation angle coincides with the angle between the resulting velocities at the relativistic addition of two velocities in the direct and reverse order.

3. The Thomas precession of the spin polarization vector at the particle motion along a curvilinear trajectory is studied. It is shown that Thomas precession is the effect of the relativistic kinematics connected with the effect of the relativistic spin rotation; the expression for the angular velocity of the Thomas precession is derived using the formula for the angle of the relativistic spin rotation.

4. The spin precession of a relativistic charged particle in the external electromagnetic field is considered. The cases of the precession in the homogeneous magnetic field and in the electric field are discussed. It is established that the angular velocity of the spin precession can be presented as a simple sum of the angular velocities of the dynamical precession due to the direct interaction of the particle magnetic moment with the magnetic field in the instantaneous rest frame and, the kinematical Thomas precession connected with the motion of a charge particle under the action of the Lorentz force.

5. The effect of the relativistic spin rotation for a system of two free particles is investigated setting their spin states in the corresponding particle rest frames. The transition from the c.m.s. of two particles to the laboratory frame is considered. It is shown that the angles of the relativistic spin rotation of two particles are generally different, except for the case of colinear vectors of particle velocities in the laboratory.

6. As a result, the effect of the relativistic spin rotation leads to the dependence of the total two-particle spin composition (the singlet and triplet fractions in the particular case of spin-1/2 particles) on the concrete frame in which the two-particle system with a nonzero vector of relative velocity is analyzed; the total spin is not a Lorentz invariant. The physical origin of this dependence is the violation of the parallelism of the spatial axes of the particle rest frames, except for the case when the Lorentz transformations to these frames are done along the directions colinear with the relative velocity (e.g., from the c.m.s. of the two particles).

This work was supported by GA Czech Republic, Grant. No. 292/01/0779, by Russian Foundation for Basic Research, Grants No. 01-02-16230 and 03-02-16210, and within the Agreements IN2P3-ASCR No. 00-16 and IN2P3-Dubna No. 00-46.


text:
# On maturation of crack patterns

## I Introduction

Consider a thin layer of a solid material elastically attached to a substrate. If the material contracts (or the substrate expands), elastic stresses appear in it. When these stresses are sufficiently high, cracks can appear in the material, giving rise to a fragmentation process. Well known examples of this phenomenon are cracking on mud and paints. In these cases the water evaporation produces the contraction of the material that is responsible for cracking. In other cases, as in the cracking of ceramic coatings, it is typically the contraction upon cooling that generates the same phenomenon. Fragmentation is known to produce a two dimensional pattern of cracks whose statistical properties have been studied theoretically, experimentally and numerically. With some variations depending on the particular case, these crack patterns are hierarchical structures, with younger cracks meeting older ones perpendicularly. Then most crack joints are ‘T’ shaped, with the horizontal part being older than the vertical part.

There is however a small number of remarkable cases in which fragmentation crack patterns undergo a ‘maturation’ process. This means that starting from a hierarchical pattern as described above, cracks can adapt smoothly to optimize its configuration. This optimization process is driven by the tendency of the crack pattern to reduce its mechanical (elastic plus crack) energy. Special conditions have to be fulfilled for this maturation to take place. To modify a given crack pattern, cracks should be able to displace laterally, and this implies typically the surmounting of enormous energy barriers (although the final state has lower energy than the original one). Particular conditions make this lateral displacement possible in (at least) two different cases.

One is the case of crack patterns formed on the ground of very cold regions of the earth, and also in other planets . In this case the frozen ground (named ‘permafrost’) cracks when the rapidly fallen temperatures of winter make the surface contract with respect to lower parts of the terrain. This first crack pattern is of the kind described above. The cracks get filled with new ice and debris, and when temperature rises after winter the cracks tend to close. However, the new material that filled the cracks is weaker than the old permafrost, and the next year cracks open almost on top of the ‘scars’ of first year cracks. However, small lateral variations can occur from one year to the next. There are many reasons that can make a crack to be shifted laterally in one direction or the other, from one year to the next. Most of these reasons (as for instance inhomogeneities in the materials) are not expected to bias the shift of the crack in one particular direction. But there is at least one reason for a crack to shift in a particular direction, and that is the tendency to reduce the energy of the crack pattern. In fact, from a statistical point of view it is reasonable to expect the crack pattern to adapt in order to reduce its energy. This tendency provides a bias for the evolution of the crack pattern in permafrost that over thousands of years is able to qualitatively modify its appearance. In fact, after maturation, crack joints become more ‘Y’ shaped, as this form has lower energy than the ‘T’ shaped original joints.

The second, better known and more remarkable example of crack pattern maturation takes place in the case of columnar jointing. It occurs in basaltic rocks when they cool after its expulsion in a volcanic event, and also in desiccating starch driven by the shrinkage due to humidity loss. In both realizations, a superficial pattern of cracks very much like the one described in the first paragraph first develops in the material. But here, this crack pattern penetrates the material as deeper parts of it cool (or dessicate). It is this progression into the interior that allows the maturation of the pattern to take place, now as a function of depth, reaching a polygonal structure whose further advance defines prismatic columns. In this case there is no true lateral movement of the cracks. But a description in terms of lateral movement can be given if we choose a reference system that moves with the penetrating crack front.

We use here a recently developed model of fracture  to describe crack patterns in a two dimensional material elastically coupled to a substrate. In the original formulation of this model cracks have to be pinned in some way in order to avoid them to move laterally (since typically this movement is unphysical). Here instead, we take advantage of this movement (driven by the tendency to minimize the energy of the system) to observe how an originally disordered pattern becomes polygonal during its maturation. We also investigate the way in which a stable polygonal pattern is modified when the degree of contraction is modified. We observe that some individual cracks disappear (terminate, in the 3D language of columnar jointing) when contraction is reduced, giving rise to local rearrangements in the pattern. This mechanism provides a way to change the mean width of the columns as a function of depth in the basalt formations and in starch, and it has been observed to occur in this last case. We finish with a discussion on what the typical width of columns in three dimensional formations is.

## II The numerical model

We use a technique recently developed  to treat fracture and cracks in the context of phase field modeling. The free energy of the system is written in terms of the strain tensor $`\epsilon _{ij}1/2(u_i/x_j+u_j/x_i)`$, with $`𝐮(𝐫)`$ being the local displacement field. We choose the form of the free energy in such a way that it reduces to the normal elastic energy for small strains, but for large strains it is able to describe cracks. This is achieved by a saturation of the free energy of the system for large values of $`\epsilon _{ij}`$. The inclusion in the free energy of terms proportional to gradients of $`\epsilon `$ produces a smoothing of cracks, which although artificial, is however very important to us. On one hand it makes the description isoptropic and insensitive to the numerical mesh we use in the calculation (as long as the discretization is much thinner than the smoothing distance of the fracture). On the other hand it allows the cracks (that in the regularized theory could be pictorially described as ‘solitons’) to move around the system to find configurations of lower energy. This wandering will model the maturation of the crack pattern. The free energy is taken to be rotationally invariant, in order to describe cracks in an isotropic material.

To write down explicitly the equations we actually solved in our two dimensional geometry, we first introduce the following notation for the independent components of $`\epsilon `$

$`e_1`$ $``$ $`(\epsilon _{11}+\epsilon _{22})/2`$ (1)
$`e_2`$ $``$ $`(\epsilon _{11}\epsilon _{22})/2`$ (2)
$`e_3`$ $``$ $`\epsilon _{12}=\epsilon _{21}`$ (3)
which are named respectively the dilation, deviatoric and shear components. These three variables are not independent. They satisfy the St. Venant compatibility constraint

$$(_x^2+_y^2)e_1(_x^2_y^2)e_22_x_ye_3=0.$$
(4)
The free energy density is

$$F(\epsilon )=\frac{F^0(\epsilon )g}{\left[1+F^0(\epsilon )/f^0\right]}$$
(5)
where

$$F^0(\epsilon )=B(e_1e_1^0)^2+\mu \left((e_2e_2^0)^2+(e_3e_3^0)^2\right)$$
(6)
and $`B`$ and $`\mu `$ are related to the two dimensional bulk and shear modulus of the material. $`e_i^0(𝐫)`$ are externally controlled functions that allow to prescribe the locally preferred state of the system, and $`g(𝐫)`$ is another (positive) function that will be used to model some random inhomogeneities in the system. The limiting value $`f_0`$ of $`F`$ for $`\epsilon \mathrm{}`$ (assuming $`g=1`$), is related to the crack energy in the model.

Regularization of cracks is provided by a gradient term $`F_g`$ in the free energy density, that we choose to be of the form

$$F_g=\underset{i=1,2,3}{}\alpha _i(e_i)^2,$$
(7)
where we have to choose $`\alpha _2=\alpha _3`$ to retain rotational invariance.

An additional ingredient that has to be added here with respect to the basic model of Ref.  is the inclusion of the elastic energy density $`F_{el}`$ of the system attached to the substrate. In terms of the displacement variables $`𝐮`$, this elastic energy can be written in the form

$$𝑑𝐫^2F_{el}\gamma 𝑑𝐫^2|𝐮(𝐫)|^2$$
(8)
where $`\gamma `$ measures the stiffness of the interaction with the substrate. As we take the components of $`\epsilon `$ to be our basic variables, we have to recast this energy in terms of them. This can be easily done in the Fourier space, and the result is

$$𝑑𝐫^2F_{el}=\gamma 𝑑𝐤^2\frac{|\stackrel{~}{e}_2(𝐤)|^2+|\stackrel{~}{e}_3(𝐤)|^2}{k^2},$$
(9)
Where $`\stackrel{~}{e}_i(𝐤)`$ are the Fourier transforms of the original $`e_i(𝐫)`$. The equations of motion are taken to be of the overdamped form, namely

$$\frac{e_i(𝐫)}{t}=\lambda \frac{\delta F}{\delta e_i(𝐫)}(i=1,2,3)$$
(10)
where

$$F=𝑑𝐫^2\left(F+F_{el}+F_g\right)$$
(11)
The Saint Venant constraint is implemented by using a Lagrange multiplier.

## III Results

We did the simulations on a square mesh of 512$`\times `$512 elements, using periodic boundary conditions. Starting from the flat configuration $`e_1(𝐫)=e_2(𝐫)=e_3(𝐫)=0`$ we simulated a uniform and abrupt contraction of the system by taking $`e_1^0(𝐫)=c`$, $`e_2^0(𝐫)=0`$, $`e_1(𝐫)=0`$, $`c=0.7`$. We introduce also a finite disorder, taking $`g(𝐫)`$ in (5) to be a random function on the lattice, uniformly distributed between 0.75 and 1.25. We keep $`B`$ and the mesh discretization $`\delta x`$ as scale-fixing parameters, and take $`\mu =0.5B`$, $`\alpha _i=0.5B\delta x^2`$ ($`i=1`$, 2, 3), $`\gamma =0.0025B/\delta x^2`$, and $`f_0=0.5B`$. The crack energy per unit length $`\eta `$ is then $`\eta =f_0\delta x`$. Under these conditions we solve the evolution equations. We see in Figs. 1 and 2 snapshots of the time evolution of the system. The figures are done by marking the points in which $`e_1>0.5`$. According to our definition of the free energy (5) and (6), this is a reasonable definition of ‘broken’ elements (results do not depends strongly on the threshold value used). In this way we are basically plotting the cracks present in the system. It is important to note that due to the finite value of $`\alpha _i`$, broken elements do not form strictly one-dimensional ‘strings’ in the system, but they clusterize, making cracks acquire a finite width, as it is apparent in the pictures. This is a crucial point to simulate an isotropic system.

We can distinguish two different stages in the temporal evolution. During the nucleation stage (Fig. 1) cracks appear rather disorderly in the system and propagate around. The pattern that forms is very dependent on many details of the simulation, as for instance the amount of disorder present. This is the kind of pattern we have described in the introduction as a fragmentation pattern. During a second stage the maturation of the pattern occurs (Fig. 2). This is observed as a progressive lateral displacement of the cracks towards a configuration of lower energy. It is necessary to emphasize again that in standard fragmentation processes this maturation cannot take place, as cracks are rigidly located in their positions. In our numerical model cracks can in fact move laterally, since this does not imply the surmounting of a large energy barrier.

The lateral movement of cracks in our model is however rather slow compared to its nucleation, and that is why it is not seen on the timescale of the nucleation stage. We stress that we are not forcing the crack pattern to become polygonal, or cracks to terminate onto other cracks, it is the system itself that prefers this kind of configuration as this reduces its energy. The final, stable pattern is that at the bottom right of Fig. 2. It corresponds to a relative minimum of the energy of the system, the absolute minimum being a perfect hexagonal pattern with a polygon size (calculated numerically with the same model) as indicated also on Fig. 2. The mature pattern contains mostly polygons of five, six, and seven sides, and a small number with four, and eight sides. They are statistically very similar to those in real columnar formations (see Fig. 4 below and Fig. 8 in Ref. ). We note that the mean area of polygons for different number of sides follows a linear relation, known as the Lewis law, after he encountered it in other two dimensional patterns. This law follows if the pattern is assumed to be maximally random.

The present results can be compared with those obtained previously using a phenomenological model for the energy of the cracked material. The present approach is however much more general than that in Ref. . Here, we are not assuming any phenomenological form of the energy as a function of the areas of the polygons, the energy of the system builds up from the free energy presented in the previous section. In addition, crack segments are not forced here to be straight, and in fact we can see in the last panel of Fig. 2 that some of them are slightly curved. The curvature occurs particularly when there is a large difference between the areas of polygons on both sides of the crack segment, always curving it in the direction in which areas tend to be closer. The reason for this is again energetic: slightly curving a crack does not pay much crack energy, but produces a gain in elastic energy if the areas of the two adjacent polygons tend to become closer to each other. This curvature has been in fact observed to occur in a full three dimensional calculation for a simple geometry.

An interesting problem to be investigated with the present model is the way in which a stable polygonal pattern changes when there are changes in the parameters that control the extent of contraction. As an outcome of this analysis we will get an idea of the expected evolution of the patterns down in the columnar formation (after the first maturation), since the thermal stresses in deeper parts of the material are lower than close to the surface. Since a lower grade of contraction corresponds to an ideal pattern with larger polygons, we may wonder what is the way (if any) in which one of our patterns adapts to the new conditions. We present in Fig. 3 the results of simulations when the extent on contraction $`c`$ is reduced. We see that there is an increase in the mean area of the polygons when $`c`$ is reduced. The area increase is not homogeneous over all polygons, but occurs due to the disappearance of particular crack segments, merging two (or three) adjacent polygons into one. After the disappearance of the crack segments there is a local rearrangement of the pattern which adjusts to the new configurations. Those regions in which no crack disappear remain perfectly stable despite the change in the contraction.

The evolution of the mean area $`A`$ of the polygons as a function of $`c`$ is shown in Fig. 4 along with the ideal area $`A^{id}`$ of the hexagons in the perfect hexagonal pattern of minimum energy. We see that the evolution tends to follow that of the ideal structure, although $`A`$ is always smaller than $`A^{id}`$ . We note also that in the present model there is a critical value of $`c`$ ($`0.42`$) for which $`A^{id}`$ diverges, and we expect the same occurs for non ideal patterns. This happens because the elastic energy per unit area gained when generating a polygonal crack pattern decays very rapidly when the size of the hexagons increases sufficiently. The sum of this elastic energy plus fracture energy may not have a minimum with respect to the area of polygons if the degree of contraction $`c`$ is too small. Note that the same does not apply to a real three dimensional columnar case (see next section).

It is remarkable that the statistical properties of the pattern do not change appreciable during this relaxation stage. In Fig. 5 we see that despite a change in the mean area by more than fifty percent, statistical distribution of polygons by number of sides and areas remain constant within numerical fluctuations associated to the finite size of the system. Note that it is precisely the ‘imperfection’ of the crack pattern that makes possible the adaptation of the mean area to a condition of lower contraction. For a perfect hexagonal pattern it is impossible to find a way to adapt the pattern slightly and obtain another hexagonal pattern with slightly larger polygon area. In our case the mean area of the pattern is increased by making some crack segments between polygons disappear.

## IV The problem of the typical column width in columnar formations

The two dimensional model we have studied is perfectly well defined, and provides values for the size of the polygons in the ideal hexagonal pattern that minimizes the total energy of the system. We want to comment at this point to what extent these two dimensional results can be applied to the full three dimensional columnar problem. For a straightforward application to be possible, the elastic energy stored in a columnarly cracked three dimensional material should be stored in a layer around the crack front of thickness $`w`$, in such a way that $`w`$ is much smaller than the typical column width $`l`$ ($`lA^{1/2}`$). If this is satisfied, the two dimensional description is directly applicable. The only consideration to be made is that two-dimensional variables have to be scaled from the three dimensional variables using $`w`$. For instance, the effective crack energy per unit length $`\eta `$ and elastic constants $`B`$ and $`\mu `$ of the two dimensional description are obtained from the real three dimensional values as $`w\eta ^{(3d)}`$, $`wB^{(3d)}`$ and $`w\mu ^{(3d)}`$. Unfortunately the condition $`w<<l`$ is never satisfied. In fact, the elastic energy of a columnar formation is stored in a portion of thickness $`wl`$ around the crack front. The coincidence of the statistical properties of our two dimensional patterns and those in true three dimensional cases indicates that these properties are robust with respect to this difference. However, the calculation of the ideal size of the perfect hexagonal pattern (and then an estimation of the typical size of non-perfect real patterns) has to be reconsidered for the three dimensional case. In fact, in our two dimensional model, in which a layer of material is attached to a substrate, there is an ideal hexagonal pattern of well defined polygon size that minimizes the energy of the system. The application of the same principle of minimizing the total energy leads in the three dimensional case to nonsense: the contribution of the fracture energy to the total energy is always much larger than the elastic contribution. A minimum can only be obtained with no cracks at all.

The correct way to pose the problem of the typical size of the polygons in three dimensions is the following: in the three dimensional case, we have a temperature profile that we assume to be dependent only on depth $`z`$, passing more or less steeply from $`T_0`$ at $`z\mathrm{}`$ to $`T_1>T_0`$ at $`z\mathrm{}`$, in such a way that a temperature front can be defined. At time $`t=0`$ the temperature front is assumed to be located at $`z=0`$. We will consider the idealized case in which the temperature profile is rigidly displaced towards the interior as a function of time, with some fixed velocity $`v`$, namely $`T(z,t)=\stackrel{~}{T}(zvt)`$. We assume that a stable polygonal pattern of fractures has formed, and that its front is located at some depth $`z_0`$, which moves down locked to the temperature profile, namely $`z_0=z_r+vt`$, where $`z_r`$ measures the relative position of the crack pattern and the temperature profile. The value of $`z_r`$ depends mainly on the typical size of the pattern $`l`$ and the overall temperature difference $`\mathrm{\Delta }T=T_1T_0`$. A previous stability analysis has shown in a simplified case that under the present conditions, patterns with different $`l`$ can propagate in a stable manner, with $`z_r`$ being a decreasing function of $`l`$, namely larger patterns are more retarded with respect to the temperature front. However there is a limit to this stable propagation. If $`l`$ or $`\mathrm{\Delta }T`$ are too small, the crack front becomes unstable: not all cracks can propagate. It is tempting to argue (and this is also based upon what is observed in three dimensional starch samples, see below) than in this case some crack segment will remain halted, and the rest of the pattern propagates. In this way $`l`$ is effectively increased and the crack front becomes closer to the temperature front, in such a way that the new pattern is now stable. In a situation in which $`\mathrm{\Delta }T`$ decays smoothly with time (whereas at the same time the front penetrates the material), we may expect that the pattern will always be located at the value $`z_{cr}`$ that marks the limit between stable and unstable propagation. This is the condition that determines the size of the columns in terms of the temperature profile, the elastic properties of the material and the crack energy. For the case of a sharp temperature jump and generalizing the two dimensional expressions for elastic and fracture energy in Ref. , we obtain that the crack front is located precisely at the border between stable and unstable regions when $`B(\alpha \mathrm{\Delta }T)^2l/\eta `$ is some constant value $`k`$ of order unity (this value is not easy to calculate). Here $`\alpha `$ is the thermal expansion coefficient and $`B`$ is a typical elastic constant of the material. From here we obtain the typical width of the columns as

$$l=k\frac{\eta }{B}(\alpha \mathrm{\Delta }T)^2$$
(12)
The typical size is then positively correlated to the crack energy, and negatively correlated with the elastic stiffness of the material, both facts being qualitatively reasonable. The size $`l`$ is also proportional to the negative second power of the temperature jump responsible for cracking. We should keep in mind however that this result is valid only for the assumed sharp step form of the temperature profile. In other cases we should search for the critical position of the crack front $`z_{cr}`$ along the lines used in Ref. . Note that $`\alpha \mathrm{\Delta }T`$ plays in three dimensions the role of the degree of contraction $`c`$ in our simulations. The preceding formula indicates that the typical size of polygons diverges only when $`\alpha \mathrm{\Delta }T0`$, contrary to the critical value of $`c`$ we found in two dimensions (see Fig. 4). This would indicate that if the driving force for cracking is slowly reduced when going deeply into the material (as may occur be due to the higher difficulty to expel heat -or humidity in starch- through the upper material) the size of the pattern should adapt by increasing their typical size, but it would never stop abruptly.

Recent tomography experiments in starch samples  show that termination and rearrangement of cracks seem to be in fact the main mechanism by which the polygonal pattern evolves in depth. In starch samples the humidity gradients are expected to be reduced when going deeper into the sample, and that is why the typical width of the columns tends to increase. However, a quantitative verification of a relation like (12) (or the equivalent one for a more realistic time dependent temperature or humidity profile) is not possible at present as it would require the in situ determination of the temperature profile under which the cracks form, and not only the observation a posteriori of column thickness as a function of depth.

## V Conclusions

We have studied numerically the formation and maturation process of a two dimensional crack pattern that is allowed to adapt to find configurations of minimum energy. The original cracks appear in a rather disordered way, but the pattern naturally evolves towards a polygonal configuration with well defined statistical properties. We argue that this maturation process occurs in cracks patterns on the ground of arctic regions (permafrost) and effectively in the columnar jointing of basalts and starches, as a function of depth. Our model allows also to study the evolution of mature polygonal patterns when the extent of contraction is reduced. We have found that in this case the pattern adapts by closing (‘terminating’ in the three dimensional interpretation) some cracks and rearranging those cracks in the immediate neighborhood. This evolution has been recently observed to occur in starch samples. Although it does not contain all features of the full three dimensional problem, our approach produces patterns of very good statistical agreement with real ones. The issue of the typical scale of the three dimensional pattern is beyond the reach of the two dimensional model, and we have provided for this case a plausible description that relates the typical size of the polygons with the elastic and thermal properties of the material, and with the details of the temperature profile.