text:
# System size dependence of strangeness production at 158 AGeV

## 1 Motivation

An increase of the ratio of strange to non-strange particles is observed in all A+A collisions compared to N+N reactions; it is found to depend on energy , size of the interacting nuclei and collision geometry . This strangeness enhancement was proposed as a signal for a deconfined state of matter. But despite rich experimental data and many theoretical efforts the origin of strangeness enhancement is still not fully understood.

In this contribution we concentrate on the dependence of strangeness enhancement on the system size and, in particular, on the collision geometry. Macroscopic, thermodynamic models relate the increase of strangeness production with system size to the transition from the canonical to the grand-canonical ensemble . The resulting disappearance of what is called canonical strangeness suppression already in rather small reaction volumes is based on a mechanism that is still open to discussion. In a microscopic picture successively excited nucleons create collision, string and energy densities varying under different experimental conditions. This provides the opportunity to search for the microscopic parameters responsible for enhanced strangeness production.

## 2 Experiment and results

NA49 is a large-acceptance hadron spectrometer . Charged pions and kaons are identified by means of energy loss and momentum information. Short-lived particles as the $`\varphi `$ or the $`\mathrm{\Lambda }`$ are measured through their hadronic decay channels, i.e. $`\varphi K^+K^{}`$ and $`\mathrm{\Lambda }p\pi ^{}`$. For the $`\mathrm{\Lambda }`$ the $`V_0`$-decay topology is used in addition. All data are corrected for acceptance, kaon decay in flight and the vertex resolution.

The data on p+p and centrality-selected Pb+Pb interactions were presented elsewhere . The light systems C+C and Si+Si were investigated using a secondary beam of Pb-fragments (for further details see ). Central collisions were chosen by setting an upper threshold on the forward energy as measured in the zero-degree calorimeter. The percentage of the inelastic cross section selected this way was used to calculate the mean number of participants $`N_{part}`$ as well as the mean number of collisions $`\nu `$ within the VENUS model (version 4.12) . In C+C (Si+Si) the 17.5% $`\pm `$ 1.5% (12.5% $`\pm `$ 1.5%) most central events correspond to 16 $`\pm `$ 1 (41.5 $`\pm `$ 1.5) participants which undergo 1.7 (2.2) collisions on average.

Transverse mass distributions in central C+C and Si+Si interactions show an exponential shape. The mass dependence of the inverse slopes (fig. 2) indicates the presence of collective transverse flow which itself increases with system size. Rapidity spectra of mesons (fig. 2) are fitted by a Gaussian to extract full yields. In the figures only statistical errors are shown; additional systematic errors for the yields amount to about 5% (10%) for pions (kaons, $`\varphi `$). The rapidity distribution of the $`\mathrm{\Lambda }`$ is rather flat over the whole measured range. A comparison to scaled p+p data shows an increase on the order of two at midrapidity which reflects both, larger stopping and strangeness enhancement.

## 3 Discussion and Interpretation

Relative strangeness production in nuclear collisions is approximately measured by the $`K/\pi `$ ratio, since roughly 70% of the produced strangeness appears in the kaons. The $`\varphi `$-meson is of interest because of its hidden strangeness. Ratios of these particles are shown in figure 3(a) as function of the number of participants. Strangeness enhancement relative to p+p interactions which increases with system size is already observed in small systems. However, the ratios are higher in central collisions of small systems than in peripheral Pb+Pb at the same number of participants indicating an important effect of the collision geometry on strangeness production. Previously  this was parametrized with the macroscopic geometrical variable $`Rb/2`$ representing the surface per volume ratio of the system. Here, this geometry effect will be discussed in terms of a microscopic reaction picture: The strategy is to search for a common scaling parameter which should give insight into the underlying reaction mechanism.

In e++e- and p+p collisions an approximate energy-independence of the $`K/\pi `$ ratio is observed in the range of interest here . This suggests that, aside from at most a few percent, higher excitation of nucleons by sequential N+N interactions in A+A collisions alone does not cause the strangeness enhancement.

Another new feature in A+A reactions compared to p+p is rescattering of secondaries. An indication that this can be neglected in light systems comes from the small collective transverse flow observed in C+C and Si+Si in comparison to central Pb+Pb collisions (fig. 2). In addition UrQMD111The author thanks the UrQMD collaboration for providing the preliminary version 1.3.  shows that the increase of the $`K/\pi `$ ratio from rescattering is minor in small systems, and, more importantly, that more rescattering takes place in peripheral Pb+Pb than in central C+C and Si+Si reactions. This is plausible because of the smaller particle multiplicity in the latter systems. Thus, if rescattering would be the source of strangeness enhancement peripheral Pb+Pb should give larger ratios.

Assuming that energy loss and rescattering are not the main sources of strangeness enhancement other features of the primary inelastic N+N collisions must play an important role. A first obvious attempt is to investigate the role of multiple interactions. As shown in , the mean number of collisions per projectile $`\nu `$ as calculated within the Glauber approach does not provide scaling; this holds also for C+C and Si+Si collisions. However, the distribution of $`\nu `$ for C+C or Si+Si and Pb+Pb reactions at the same $`K/\pi `$ ratio shows that the fraction of nucleons which undergo multiple collisions is similar. In fact scaling is observed in this latter variable (fig. 3(b)).

The sequential N+N reactions happen in close vicinity to each other, thus causing a high density of interactions in space and time. Indeed, the mean space-time density of all inelastic collisions during the first phase $`\rho _{\text{inel coll}}`$, i.e. when the nucleons pass through each other, as calculated in the center of mass system of the collision within the UrQMD model (version 1.3) serves also as scaling parameter (fig. 3(c)).

## 4 Conclusions

The presented data show that the number of participating nucleons in A+A collisions is not the decisive variable for strangeness enhancement (fig. 3(a)). However, taking advantage of this observation it is found that features connected to the primary inelastic N+N collisions provide a common description of the $`K^+/\pi ^\pm `$, $`K^{}/\pi ^\pm `$ and the $`\varphi /\pi ^\pm `$ ratio in various A+A collision systems. This implies that sequential N+N interactions within a small volume and period of time are not independent of each other with regard to strangeness production.

The collision density as defined here is on the one hand related to the energy density in the system. On the other hand a high density of interactions may also mean overlap of strings thus leading to string fusion and enhanced string tension. Both, string overlap and energy density have often been discussed as being relevant for enhanced strangeness production. In addition a volume in which a high energy density or overlapping strings and thus a high space-time density of successive collisions exist might decay in a quantum-mechanically coherent fashion, and thus represent the hadronizing volume implied by the canonical and grand-canonical versions of the statistical hadronization model.

Acknowledgements: This work was supported by the Director, Office of Energy Research, Division of Nuclear Physics of the Office of High Energy and Nuclear Physics of the US Department of Energy (DE-ACO3-76SFOOO98 and DE-FG02-91ER40609), the US National Science Foundation, the Bundesministerium für Bildung und Forschung, Germany, the Alexander von Humboldt Foundation, the UK Engineering and Physical Sciences Research Council, the Polish State Committee for Scientific Research (2 P03B 130 23 and 2 P03B 02418), the Hungarian Scientific Research Foundation (T14920 and T32293), Hungarian National Science Foundation, OTKA, (F034707), the EC Marie Curie Foundation, the Polish-German Foundation, and Bergen Computational Physics Laboratory in the framework of the European Community - Access to Research Infrastructure action of the Improving Human Potential Programme.


format: text

text:
# One-Dimensional Confinement and Enhanced Jahn-Teller Instability in LaVO3

\[

## Abstract

Ordering and quantum fluctuations of orbital degrees of freedom are studied theoretically for LaVO3 in spin-C-type antiferromagnetic state. The effective Hamiltonian for the orbital pseudospin shows strong one-dimensional anisotropy due to the negative interference among various exchange processes. This significantly enhances the instability toward lattice distortions for the realistic estimate of the Jahn-Teller coupling by first-principle LDA+$`U`$ calculations, instead of favoring the orbital singlet formation. This explains well the experimental results on the anisotropic optical spectra as well as the proximity of the two transition temperatures for spin and orbital orderings.

\]

Orbital degrees of freedom are playing key roles in magnetic and charge transport properties of transition metal oxides . Especially it has been recognized that, in these strongly-correlated systems, spatial shapes of orbitals can give rise to an anisotropic electronic state even in the three-dimensional (3D) perovskite structure . There the spin ordering (SO) and orbital ordering (OO) are determined self-consistently .

Perovskite vanadium oxides, AVO3 (A is rare-earth element), are typical $`t_{2g}`$ electron systems which show this interplay between orbital and spin degrees of freedom . Both magnetic and orbital transition temperatures, $`T_\mathrm{N}`$ and $`T_\mathrm{o}`$, respectively, change systematically according to the ionic radius of the A atom , which controls the bandwidth through the tilting of VO6 octahedra. For smaller ionic radii (smaller bandwidth) such as A=Y, $`T_\mathrm{o}`$ for OO of G-type (3D staggered) is much higher than $`T_\mathrm{N}`$ for SO of C-type (rod-type) . As the ionic radius increases (the bandwidth increases), $`T_\mathrm{o}`$ decreases while $`T_\mathrm{N}`$ increases, and finally they cross between A=Pr and Ce . In LaVO3, the SO occurs at $`T_\mathrm{N}143`$K first, and at a few degrees below $`T_\mathrm{N}`$ the OO takes place . A remarkable aspect here is its proximity of $`T_\mathrm{N}`$ and $`T_\mathrm{o}`$, which is also observed for all the compounds with $`T_\mathrm{N}>T_\mathrm{o}`$, i.e., CeVO3  and La1-xSrxVO3 ($`x<0.17`$) . Therefore, in LaVO3, the magnetic correlation appears to develop primarily and to induce the orbital transition immediately once the SO sets in.

Another interesting aspect of LaVO3 is the large anisotropy in the electronic state, which has recently been explored by the optical spectra . Figure 1 shows the temperature dependence of the spectral weights, $`I_c`$ along the $`c`$ direction and $`I_{ab}`$ within the $`ab`$ plane, which is obtained from the data in ref. REFERENCES. Here we define the spectral weight as an integration of the optical conductivity up to the isosbestic (equal-absorption) point at $`2.8`$eV, namely, $`I_\mu \frac{2m_0}{\pi e^2n}_0^{2.8\mathrm{eV}}\sigma _\mu (\omega )𝑑\omega `$ ($`\mu =c`$ or $`ab`$), where $`m_0`$ and $`n`$ are the free electron mass and the density of V atoms, respectively. The most striking feature is the temperature dependence. $`I_c`$ grows rapidly below $`T_\mathrm{N}`$ while $`I_{ab}`$ is almost temperature independent. Therefore the temperature dependence is almost 1D although the ratio $`I_c/I_{ab}2`$ is not so large.

Compared with $`e_g`$ systems, the Jahn-Teller (JT) coupling in $`t_{2g}`$ systems is expected to be weak, and the quantum fluctuation and/or the singlet formation of the orbital degrees of freedom is a keen issue. Recently, Khaliullin et al. claimed that orbital singlet correlation along the $`c`$ axis is the driving force to realize the ferromagnetic spin exchange in the C-type antiferromagnetic (AF) phase in LaVO3 . On the other hand, the C-type SO state with the G-type OO has been obtained by the mean-field theory  and the first-principle calculation  which are justified for rather weakly-correlated cases and do not take account of quantum fluctuations seriously. There, the AF interactions within the $`ab`$ plane concomitant with a single occupation of $`xy`$ orbitals play a key role to stabilize the SO and OO. Hence, these two pictures are quite different. Since there are competing interactions with the orbital quantum nature, such as JT coupling and 3D orbital exchange couplings, it is highly nontrivial to what extent the quantum fluctuations are important under the realistic situations, nevertheless the quantitative study on this issue has been missing.

In this Letter, we study ordering and fluctuations of orbital degrees of freedom in the C-type AF phase in LaVO3. An effective orbital model including the JT and the relativistic spin-orbit couplings is derived using the parameters obtained by first-principle LDA+$`U`$ calculations and the optical experiments. We show that the model exhibits a strong 1D anisotropy which explains well the experimental results for the optical spectra and the proximity of $`T_\mathrm{o}`$ to $`T_\mathrm{N}`$. It is concluded that the enhanced JT instability due to the 1D confinement dominates the orbital singlet formation in LaVO3.

Now we derive the effective orbital model. We start from the strong-coupling limit of the Hubbard model with three-fold orbital degeneracy for the $`t_{2g}`$ orbitals . The system contains two $`d`$ electrons at each V atom which form the high-spin $`S=1`$ state due to the Hund’s-rule coupling. To focus on the orbital sector in this spin-orbital coupled Hamiltonian, we assume the C-type SO. This is a reasonable approximation because the spin $`S=1`$ has less quantum nature compared to $`S=1/2`$ and the C-type SO is obtained by the mean-field calculation  and the first-principle calculation  as mentioned above. At the same time, we assume that the $`xy`$ orbital is singly occupied at each V atom, which drives the AF coupling in the $`ab`$ plane . Then the second electron goes to either $`yz`$ or $`zx`$ orbital. We assign a pseudospin state $`\tau ^z=\pm 1/2`$ for the occupancy of these two orbitals . Finally, our total Hamiltonian is written as $`=_{\mathrm{orb}}^c+_{\mathrm{orb}}^{ab}+_{\mathrm{JT}}+_{LS},`$ each term of which is described below.

The orbital exchange term for the $`c`$ direction is

$$_{\mathrm{orb}}^c=\underset{ij}{}J_c(\stackrel{}{\tau }_i\stackrel{}{\tau }_j1/4),$$
(1)
while that within the $`ab`$ plane is

$$_{\mathrm{orb}}^{ab}=\underset{ij}{}\left(J_{ab}^{}\tau _i^z\tau _j^zJ_{ab}^+/4J_{ab}^{xy}/4\right).$$
(2)
Here the summations are taken for nearest-neighbor pairs. The coupling constants are given by

$`J_c`$ $`=`$ $`4t^2/(U^{}J_\mathrm{H}),`$ (3)
$`J_{ab}^\pm `$ $`=`$ $`{\displaystyle \frac{2t^2}{3(U^{}J_\mathrm{H})}}+{\displaystyle \frac{4t^2}{3(U^{}+2J_\mathrm{H})}}\pm {\displaystyle \frac{t^2}{U+2J_\mathrm{H}}}\pm {\displaystyle \frac{t^2}{U}},`$ (4)
$`J_{ab}^{xy}`$ $`=`$ $`4t^2/(U+2J_\mathrm{H})+4t^2/U,`$ (5)
where $`U,U^{}`$ and $`J_\mathrm{H}`$ are the intra-orbital, the inter-orbital Coulomb interaction and the Hund’s-rule coupling, respectively. Neglecting the small tilting of VO6 octahedra in LaVO3, the transfer integral $`t`$ is taken to be diagonal which strongly depends on the direction, for instance, $`t_{ij}^{yz}=t_{ij}^{zx}=t`$ and otherwise zero in the $`c`$ direction. From the analysis in ref. REFERENCES, the parameters are estimated as $`U2.25`$eV, $`U^{}1.93`$eV and $`J_\mathrm{H}0.16`$eV. The transfer integral $`t`$ is set to be $`0.12`$eV based on the estimate of the bandwidth $`1`$eV in first-principle calculations . Then the orbital exchange interaction along the $`c`$ axis is estimated as $`J_c33`$meV, while $`J_{ab}^{}2`$meV. We set $`J_c=1`$ as an energy unit in the following calculations.

Here we point out two important features in eqs. (1) and (2). One is that the exchange in the $`c`$ direction is Heisenberg-type while that within the $`ab`$ plane is Ising-type. From this, one might expect a strong quantum fluctuation in the $`c`$ direction as pointed out in ref. REFERENCES. However, this quantum nature becomes relevant only when the JT coupling is negligibly small, and this is not the case in LaVO3 as discussed in the following. The other important feature is the large 1D anisotropy in the orbital exchange couplings. Note that the negative interference among different perturbation processes occurs in the in-plane coupling $`J_{ab}^{}`$ in eq. (4), which results in the ratio of the exchange couplings $`J_c/J_{ab}^{}17`$.

The JT coupling in the subspace of $`\tau `$ is given by

$$_{\mathrm{JT}}=\underset{i}{}gQ_i\tau _i^z+\frac{1}{2}\underset{i}{}Q_i^2,$$
(6)
where $`Q_i`$ is the JT phonon coordinate at site $`i`$. We neglect the kinetic energy of phonons, namely, regard $`Q_i`$ as a classical variable. It is crucial to estimate the coupling constant $`g`$, and we have done the following first-principle calculation . Assuming the tetragonal symmetry, we calculate the total energy as a function of the JT distortion as shown in Fig. 2. This gives the JT stabilization energy $`27`$meV, which approximately corresponds to $`E_{\mathrm{JT}}g^2/8`$ in our model. Thus we obtain the estimate $`E_{\mathrm{JT}}0.8J_c`$ ($`g2.6`$), which is appreciable and cannot be neglected as in ref. REFERENCES.

The last term is the relativistic spin-orbit coupling, which may be important in the $`t_{2g}`$ systems. We obtain the effective Hamiltonian by projecting the original form $`_{LS}=_i\lambda \stackrel{}{L}_i\stackrel{}{S}_i`$ to the subspace of $`\tau `$ by using the experimental fact that the spins lie within the $`ab`$ plane . Here $`\stackrel{}{L}_i`$ is the orbital angular momentum. Since $`L_x`$ has matrix elements between $`xy`$ and $`(yz,zx)`$ in this case, the spin-orbit interaction is represented by $`_{LS}=_i\zeta \tau _i^z,`$ where $`\zeta \lambda ^2/\mathrm{\Delta }`$ and $`\mathrm{\Delta }`$ is the energy separation between $`xy`$ and ($`yz,zx`$) orbital levels. This indicates that the spin-orbit coupling corresponds to the pseudo magnetic field along the $`z`$ direction. Using $`\lambda 20`$meV in V atom and the estimate of $`\mathrm{\Delta }1`$eV in the band calculation , we estimate $`\zeta 0.4`$meV$`0.01J_c`$. This is small enough to be neglected in the following calculations.

First we discuss the anisotropic temperature dependence in the spectral weight. The spectral weight is generally given by the kinetic energy of the system in the ground state $`K`$ as $`I=\pi K/4N`$ where $`N`$ is the system size and we set $`e^2=1`$ . In the strong-coupling limit, this is calculated by the exchange energy as shown below. Let us consider the Heisenberg model for the strong-coupling limit of the Hubbard model at half filling. The kinetic energy in the Hubbard model is calculated as $`K=t/t=/\mathrm{ln}t`$. Here the bracket denotes the expectation value in the ground state. Since the exchange coupling $`J`$ in the Heisenberg model is proportional to $`t^2`$, we obtain $`I=(\pi /4N)\times 2/\mathrm{ln}J=\pi _{\mathrm{Heis}}/2N`$ where $`_{\mathrm{Heis}}`$ is the Heisenberg Hamiltonian. Using this formula, we calculate the spectral weight for the present model by the exact-diagonalization in the ground state for $`4\times 4`$-site lattice embedded in the $`ac`$ plane. In Fig. 3 (a), the obtained values of $`\mathrm{\Delta }I_c\pi J_c\stackrel{}{\tau }_i\stackrel{}{\tau }_j/2N`$ and $`\mathrm{\Delta }I_{ab}\pi J_{ab}^{}\tau _i^z\tau _j^z/2N`$ are plotted as a function of the JT coupling $`g`$. Note that these exchange correlations correspond to the enhancements of the spectral weights from the high-temperature limit to the ground state. The results indicate that, the contribution in the $`c`$ axis is much larger than that in the $`ab`$ plane. The ratio becomes larger than $`20`$ for the realistic value of $`g2.6`$. This explains well the anisotropic temperature dependence of the spectral weight in Fig. 1.

We also discuss the magnitude of the total spectral weight in the ground state. Figure 3 (b) shows the calculated spectral weights $`I_\mu =\pi _{\mathrm{orb}}^\mu /2N`$ which include the contributions from the constants in eqs. (1) and (2). As shown in Fig. 3 (b), the ratio $`I_c/I_{ab}`$ is almost $`1`$ for the realistic value of $`g`$, which is much smaller than that of the temperature dependent part in Fig. 3 (a). This ratio is comparable but smaller than the experimental result $`I_c/I_{ab}2`$ in Fig. 1. This might be due to the facts that (i) the cut-off energy $`2.8`$eV for the integration is not large enough to take all the contribution of $`I_{ab}`$, and (ii) the transfer integrals for the different orbitals are different because they depend on the energy difference $`\mathrm{\Delta }_{dp}`$ between the $`d`$ and the oxygen $`p`$ orbitals as the energy denominator. Since the huge anisotropy in temperature dependences is reproduced with the moderate anisotropy in the total weights, we believe that the orbital 1D confinement in our model plays a major role in the anisotropic electronic state in this material.

We now turn to the discussion on the transition temperatures. In order to obtain the phase diagram at finite temperatures, we apply the transfer matrix method  combined with the unrestricted Hartree-Fock approximation for the interchain coupling $`J_{ab}^{}`$ . We start from a configuration of $`\{Q_i\}`$ and calculate the expectation value $`\tau _i^z`$ numerically. Then, from the self-consistent equation $`Q_i=g\tau _i^z`$ which is obtained as the energy minimization for $`Q_i`$, we define a new configuration of $`\{Q_i\}`$. We repeat this procedure until $`\{Q_i\}`$ is optimized and the energy is minimized. In this treatment, $`J_{ab}^{}`$ is effectively absorbed into the JT coupling $`g`$ as $`_i(gQ_i2J_{ab}^{}\tau _i^z)\tau _i^z=_i\overline{g}^2\tau _i^z\tau _i^z`$ where $`\overline{g}^2g^2+2J_{ab}^{}`$. The obtained transition temperature for the orbital/lattice ordering transition is plotted as a function of $`\overline{g}`$ in Fig. 4. We note that for large $`g`$, the transition temperature diverges as $`T_\mathrm{o}g^2/4`$ which corresponds to the JT energy gap. However, this is an artifact of the mean-field-type treatment, and in reality $`T_\mathrm{o}`$ should stay at a constant of the order of $`J_c`$ since the coupling by $`J_c`$ between the neighboring sites determines $`T_\mathrm{o}`$ in the limit of large $`g`$.

Considering the realistic value of $`J_c33`$meV and $`E_{\mathrm{JT}}27`$meV ($`g\overline{g}2.6`$) for LaVO3, the orbital transition temperature is estimated as $`T_\mathrm{o}800`$K which is much higher than the observed $`T_\mathrm{N}143`$K as indicated in Fig. 4. One might think that this contradicts with the experimental fact $`T_\mathrm{N}>T_\mathrm{o}^{\mathrm{exp}}`$. However the phase diagram is obtained assuming the C-type SO, which induces the 1D confinement of the orbital degrees of freedom with the enhanced $`J_c`$. Note that the disorder of the spins should reduce the effective orbital exchange as easily shown in the spin-orbital coupled Hamiltonian. It is also well-known that 1D systems have an enhanced instability to lattice distortions compared with higher dimensions. Therefore, $`T_\mathrm{o}^{1\mathrm{D}}=T_\mathrm{o}`$ under this 1D orbital confinement can be higher than $`T_\mathrm{o}^{3\mathrm{D}}`$ without the SO when the JT coupling governs the OO transitions. In real materials with $`T_\mathrm{o}>T_\mathrm{N}`$, $`T_\mathrm{o}^{3\mathrm{D}}`$ decreases as the bandwidth increases, which indicates the relevance of the JT coupling. (If the 3D orbital exchange couplings dominate, $`T_\mathrm{o}^{3\mathrm{D}}`$ should increase as $`T_\mathrm{N}`$ does.) Then, when the inequality $`T_\mathrm{o}^{3\mathrm{D}}<T_\mathrm{N}<T_\mathrm{o}^{1\mathrm{D}}`$ is satisfied, the OO transition with the JT lattice distortion should take place as soon as the SO grows and induces the 1D confinement in the orbital channel. In this scenario, comparing $`T_\mathrm{o}^{1\mathrm{D}}`$ with $`T_\mathrm{N}`$, we can estimate the lower bound for the value of $`g`$ to realize this proximity of the transition temperatures. Our estimation for this lower bound from Fig. 4 is $`g1`$ which is consistent with the estimate in Fig. 2.

Let us discuss this proximity of $`T_\mathrm{N}`$ and $`T_\mathrm{o}^{\mathrm{exp}}`$ by the Ginzburg-Landau type argument. In this spin-orbital coupled system, we have the term in which the SO parameter $`M`$ and the OO one $`O`$ are coupled as $`\{a(T)bM^2\}O^2`$. Here $`a(T)`$ is the coefficient for the second-order term for OO without SO which is given by $`a(T)=a^{}(TT_\mathrm{o}^{3\mathrm{D}})`$. Our results indicate the relation $`a(T)bM_{\mathrm{sat}}^2=a^{}(TT_\mathrm{o}^{1\mathrm{D}}),`$ where $`M_{\mathrm{sat}}`$ is the saturated magnetic moment. $`T_\mathrm{o}^{\mathrm{exp}}`$ is given by solving the equation $`a(T_\mathrm{o}^{\mathrm{exp}})bM^2(T_\mathrm{o}^{\mathrm{exp}})=0`$. Assuming $`M(T)=M_{\mathrm{sat}}\sqrt{(T_NT)/T_N}`$ for simplicity, the difference between $`T_\mathrm{N}`$ and $`T_\mathrm{o}^{\mathrm{exp}}`$ is given by $`\delta T/T_\mathrm{N}(T_\mathrm{N}T_\mathrm{o}^{\mathrm{exp}})/T_\mathrm{N}=(T_\mathrm{N}T_\mathrm{o}^{3\mathrm{D}})/(T_\mathrm{N}+\mathrm{\Delta }T),`$ where $`\mathrm{\Delta }T=T_\mathrm{o}^{1\mathrm{D}}T_\mathrm{o}^{3\mathrm{D}}`$. Considering the systematic changes of $`T_\mathrm{N}`$ and $`T_\mathrm{o}`$ for A-site ions , we expect that $`T_\mathrm{o}^{3\mathrm{D}}`$ is slightly lower than $`T_\mathrm{N}`$ in LaVO3 and CeVO3. Assuming $`T_\mathrm{o}^{3\mathrm{D}}=0.8T_\mathrm{N}`$, we plot the expected $`T_\mathrm{o}^{\mathrm{exp}}`$ in Fig. 4 as the gray curve. For the realistic value of $`g`$, we have $`T_\mathrm{o}^{\mathrm{exp}}`$ quite close to $`T_\mathrm{N}`$ as observed in these compounds.

To summarize, we have investigated the role of orbitals to understand the electronic state in LaVO3. We have derived the effective orbital model with strong one-dimensional anisotropy assuming the C-type spin ordering. We conclude that with the realistic Jahn-Teller coupling, the orbital 1D confinement leads to an enhanced instability toward lattice distortions suppressing the orbital quantum nature. This gives a comprehensive description of the anisotropy in the optical spectra and the proximity of the critical temperatures of magnetic and orbital transitions.

The authors appreciate Y. Tokura, S. Miyasaka, T. Arima, G. Khaliullin and B. Keimer for fruitful discussions. Part of numerical computations are done using TITPACK ver.2 developed by H. Nishimori, who is also thanked.


text:
# A master equation approach to option pricing

## 1 Introduction

As is well-known the seminal Black-Scholes analysis  which leads to the fair value of an option is based mainly upon the following assumptions. First the stock price, i.e. the underlying, follows a geometric Brownian motion. Second a hedge position is formed with a portfolio of short underlying and a long position of a number of European options. Then an arbitrage argument leads to the renowned Black–Scholes partial differential equation determining the value of the option. Of course, for simple cases, e.g. constant interest rate and volatility, explicit analytical solutions to the equation are known . However, for more involved cases one has to rely upon numerical methods. Here we are not concerned with deterministic methods, e.g. finite differences, but rather with Monte Carlo methods .
The idea behind the canonical Monte Carlo method is to exploit the fact that the fair value of an option is given by the present value of the expected payoff at expiry under the risk neutral measure. Thus the standard Monte Carlo approach is based upon the simulation of a geometric Brownian motion for the underlying asset until expiry. Then the payoff is computed and discounted up to the current time. By averaging over different realizations of this stochastic process the current option price is estimated . In spirit this approach is similar to the standard Monte Carlo random walk technique for the solution of partial differential equations with boundary value problems .
Monte Carlo methods offer easy to understand numerical algorithms which can be easily applied to quite complicated, e.g. path dependent or correlated multi-asset, options . In addition they allow for a straightforward inclusion of stochastic terms such as the interest rate or volatility. From a numerical point of view they are especially suited to problems with many degrees of freedom and the algorithm can be easily run in parallel.
However Monte Carlo simulations of the Black-Scholes equation are usually slower than comparable finite difference solutions of the partial differential equation. Another disadvantage of the standard Monte Carlo approach is that while the application to European options is straightforward, the valuation of American options is more involved. Using a generalisation of the canonical Monte Carlo method one has to assure that the early exercise condition stating that the option price is always above the current payoff is not violated sometime in the future. If this happens and the option value is below the payoff, the option is exercised and the option value is given by the payoff function. A Monte Carlo simulation based on the dynamics of the underlying asset has to keep track of all these points in asset time space which makes the algorithm ineffective. Nevertheless more advanced Monte Carlo algorithms are available .
Since the Monte Carlo methods described above are based on a model of the underlying dynamics of asset values these approaches could be named microscopic. In contrast now we are going to follow a different strategy, which we would like to name mesoscopic, which has already been applied with success to simulations of turbulence in fluid dynamics , to the investigation of hydrodynamic fluctuations  as well as to magnetohydrodynamics . The same approach has been shown to lead to fast Monte Carlo algorithms for the balance equations of nonequilibrium thermodynamics , for the simulation of chemical reactions  and reaction–diffusion processes .
The basic idea of the mesoscopic approach  is to regard the value of the option, say $`V`$ as the expectation value of an underlying fluctuating (mesoscopic) stochastic option value, say $`\theta `$. A master equation for $`\theta `$ is easily constructed in such a way that the expectation value of the multivariate stochastic variable $`\theta `$ satisfies the macroscopic Black–Scholes equation. The stochastic process defined in terms of a master equation is easily simulated and allows for a Monte Carlo algorithm for the direct simulation of the Black–Scholes equation.
This paper is structured as follows. In section 2 the mesoscopic master equation approach is motivated. The dynamics of the underlying asset is modelled as a piecewise deterministic process (PDP). The derivation of the fair option price for a mesoscopic stochastic option value, which parallels the Black–Scholes one, shows that the expectation value of the stochastic option price satisfies the “macroscopic” Black–Scholes equation. In section 3 a master equation for the mesoscopic fair option value is constructed, and a stochastic simulation algorithm is applied to the valuation of European and American options. In the last section perspectives and further development of the approach are indicated.

## 2 Black–Scholes Equation

The standard textbook derivation of the Black–Scholes equation assumes a geometric Brownian motion for the underlying asset which is described by the stochastic differential equation

$$ds=\mu sdt+\sigma sdW.$$
(1)
Here $`s`$ is the asset value, $`\mu s`$ and $`\sigma ^2s^2`$ denote drift and variance of the random walk. Since $`dW`$ is the increment of a Wiener process the resulting trajectory will be continuous. For more realistic models instantaneous jumps of the underlying, described by additional stochastic terms based on Poisson processes , may be included in (1). These so called jump diffusion models, are used e.g. to incorporate the effect of information about stocks which arrives at random times .

### 2.1 Black–Scholes equation from a piecewise deterministic process (PDP)

In contrast to the standard derivation of the Black–Scholes equation based on (1), this section demonstrates how the derivation can be paralleled by modelling the underlying in terms of a piecewise deterministic process . This is possible since a diffusion process can be represented as the continuous limit of an appropriate jump process. This approach, despite describing the same dynamics as the standard approach, has the advantage of clearly demonstrating that it is quite obvious to interpret the option value as a stochastic variable and thus to interpret the Black–Scholes equation as the macroscopic expectation value of a mesoscopic stochastic option value. In addition the formulation in terms of PDPs simplifies the inclusion of additional jump processes since one has to deal only with one type of stochastic process. Another important point is that for such processes one can easily proceed to a master equation formulation, for which powerful numerical algorithms exist, as will be shown later on.
Consider now the stochastic differential equation

$$ds=s\mu dt+\delta sdN^+\delta sdN^{},$$
(2)
with $`dN^\pm =s^2\sigma ^2/(2\delta s^2)dt`$. The Poisson increments $`dN^\pm `$ are either $`0`$, in which case the deterministic evolution takes place, or $`1`$, which describes an instantaneous jump of size $`\pm \delta s`$. Stochastic process of this general form are thus also known as piecewise deterministic processes .
In the limit $`\delta s0`$ this stochastic differential equation corresponds to a geometric Brownian motion with drift $`s\mu `$ and variance $`s^2\sigma ^2`$ . Hence in this limit the above equation generates the same dynamics as Eq. (1), but is formulated as a PDP.
Using the above stochastic differential Eq. (2) instead of Eq. (1) as starting point for the classical Black–Scholes derivation  one obtains in the limit $`\delta s0`$ for the expectation value of the stochastic option price $`v`$

$`{\displaystyle \frac{v(s,t)}{t}}={\displaystyle \frac{s^2\sigma ^2}{2}}{\displaystyle \frac{^2v(s,t)}{s^2}}rs{\displaystyle \frac{v(s,t)}{s}}+rv(s,t),`$ (3)
where $`r`$ denotes the riskfree interest rate. If the expectation value $`V(S,t)=v(s,t)`$ is identified with the macroscopic option price then one obtains the Black–Scholes equation for a European option on a non dividend paying asset.
Hence the traditional Black-Scholes equation can be interpreted as the expectation value of a stochastic process. Additional not infinitesimal jump processes can now easily be included, but this of course will change the hedging strategy  and the expectation value of the option price will then follow equations similar to the Black–Scholes equation.

### 2.2 Macroscopic Black–Scholes equation

The above derivation assumed a non-dividend paying stock, this assumption can be easily omitted to obtain

$`{\displaystyle \frac{V(S,t)}{t}}={\displaystyle \frac{\sigma ^2S^2}{2}}{\displaystyle \frac{^2V(S,t)}{S^2}}(rq)S{\displaystyle \frac{V(S,t)}{S}}+rV(S,t).`$ (4)
Here $`\sigma `$ and $`r`$ again denote the volatility and the interest rate and the continuous dividend is given by $`q`$. In order to numerically solve this partial differential equation the relevant variables are made dimensionless and the time direction is reversed, i.e. we have,

$$\stackrel{~}{S}=\mathrm{ln}(S/E),\stackrel{~}{V}(\stackrel{~}{S},\tau )=V(S,t)/E,\tau =\frac{\sigma ^2}{2}(Tt).$$
(5)
In this way one obtains

$`{\displaystyle \frac{\stackrel{~}{V}(\stackrel{~}{S},\tau )}{\tau }}={\displaystyle \frac{^2\stackrel{~}{V}(\stackrel{~}{S},\tau )}{\stackrel{~}{S}^2}}+(k_q1){\displaystyle \frac{\stackrel{~}{V}(\stackrel{~}{S},\tau )}{\stackrel{~}{S}}}k_r\stackrel{~}{V}(\stackrel{~}{S},\tau ),`$ (6)
where $`k_q=2(rq)/\sigma ^2`$ and $`k_r=2r/\sigma ^2`$ have been introduced. If the parameters of the Black–Scholes equation (4) are assumed to be constant in time, then the analytic solution to the above equation is available, see e.g. .

## 3 Master Equation formulation of the Black-Scholes Equation

As seen in section 2.1 the Black–Scholes equation can be interpreted as a deterministic equation which governs the expectation value of a stochastic option price. In this section we will describe how the stochastic dynamics of the option price can be formulated in terms of a master equation. For the sake of simplicity we will consider the Black-Scholes equation for a European call in the dimensionless form of Eq. (6). The generalisation to a put is straightforward and the application to American options is discussed in section 3.6.

### 3.1 General Theory

A general approach in statistical physics is to derive macroscopic equations from a known microscopic dynamics. Here the opposite approach is followed: given a partial differential equation for a macroscopic observable we construct a mesoscopic stochastic process such that the expectation value of the stochastic variable is governed by the original partial differential equation. In this section the general theory needed for the construction of a stochastic process underlying a partial differential equation is presented, for a review see .
In order to define the underlying stochastic process the state space of the system has to be given. To this end the space variable $`x`$ is discretized such that $`\theta _\lambda `$ denotes the state of the system at the discrete position $`x_\lambda ,\lambda =1,\mathrm{},n`$. Hence the state space $`\mathrm{\Gamma }`$ of the system is given by $`\mathrm{\Gamma }=\{\theta \theta ^n\}`$. The stochastic process is defined by the joint probability distribution $`P=P(\left\{\theta _\lambda \right\},t)`$ giving the probability of finding the values $`\{\theta _\lambda \}`$ at time $`t`$. The time development of the probability density $`P`$ is given by a master equation of the general form

$$\frac{}{t}P(\theta ,t)=𝒜P(\theta ).$$
(7)
The operator $`𝒜`$ is defined in terms of diffeomorphic maps $`b`$ acting on the state space $`b:\mathrm{\Gamma }\mathrm{\Gamma }`$, and corresponding operators $`_b`$ acting on functions $`F:\mathrm{\Gamma }`$ in the following way

$$_bF(\theta )=F(b^1(\theta )).$$
(8)
The most general operator $`𝒜`$ needed in this paper is of the form

$$𝒜=c\underset{\mu }{}[\text{Det}(b_\mu )]^1_{b_\mu }𝕀,$$
(9)
where $`\text{Det}(b_\mu )`$ denotes the determinant of the Jacobian matrix of the map $`b_\mu `$ and $`𝕀`$ denotes the identity. In order to generate a valid stochastic process the operator $`𝒜`$ has to fulfil (i) all transition probabilities $`w`$ are positive and (ii) $`𝒜F=0`$ in order to preserve the normalisation of the probability density.
With these definitions the macroscopic equation of motion for a general observable $`F`$ can now easily be computed, one obtains

$$\frac{}{t}F=c\underset{\mu }{}F(b_\mu (\theta ))F(\theta ).$$
(10)
The above equation enables the computation of the macroscopic expectation value of a general observable $`F`$ given an arbitrary multivariate Markov process defined by the time evolution operator $`𝒜`$.
For a master equation formulation of the stochastic process one needs the transition probability

$$w(\theta ,\stackrel{~}{\theta })=c\underset{\mu }{}\delta (\theta ,b_\mu (\stackrel{~}{\theta })).$$
(11)
Hence the canonical form of the master equation is obtained

$$\frac{}{t}P(\theta ,t)=D\stackrel{~}{\theta }\left\{w(\theta ,\stackrel{~}{\theta })P(\stackrel{~}{\theta },t)w(\stackrel{~}{\theta },\theta )P(\theta ,t)\right\}.$$
(12)
Now the approach sketched above is applied to the Black-Scholes equation. Hence the interesting price range of the underlying $`x=\stackrel{~}{S}`$ is divided in $`n`$ discrete points $`x_\lambda ,\lambda =1,\mathrm{},n`$ with distances $`\mathrm{\Delta }_\lambda =x_\lambda x_{\lambda 1},\lambda =2,\mathrm{},n`$.

### 3.2 Stochastic process with uniform discretization

The easiest possibility to construct the stochastic process is to use a uniform discretization of the state space $`\mathrm{\Delta }_\lambda =\delta l`$. The Black-Scholes equation (6) consists of a diffusive, a convective and a part corresponding to a chemical reaction. The stochastic process underlying the Black-Scholes equation will also consist of three parts corresponding to these terms in the partial differential equation. The total time evolution operator $`𝒜`$ is thus defined by

$$𝒜_{\text{BS}}=𝒜_{\text{diff}}+𝒜_{\text{conv}}+𝒜_{\text{chem}}.$$
(13)
The above time evolution operators are defined through their corresponding maps, see Eq. (9). The map

$$a_\mu ^\pm =\{\begin{array}{cc}\theta _\mu \hfill & \theta _\mu \alpha \theta _\mu \hfill \\ \theta _{\mu \pm 1}\hfill & \theta _{\mu \pm 1}+\alpha \theta _\mu \hfill \end{array},$$
(14)
describes the diffusive part, while the maps $`b_\mu `$ and $`c_\mu `$

$$b_\mu =\{\begin{array}{cc}\theta _\mu \hfill & \theta _\mu \alpha \theta _\mu \hfill \\ \theta _{\mu 1}\hfill & \theta _{\mu 1}+\alpha \theta _\mu \hfill \end{array},c_\mu =\{\begin{array}{cc}\theta _\mu \theta _\mu \alpha \theta _\mu \hfill & \end{array},$$
(15)
correspond to the convective part and to the chemical reaction term. Hence the time evolution operators for the three parts of the stochastic process are given by

$`𝒜_{\text{diff}}`$ $`=`$ $`{\displaystyle \frac{1}{\alpha \delta l^2}}{\displaystyle \underset{\mu }{}}{\displaystyle \frac{1}{1\alpha }}\{𝒜_\mu ^++𝒜_\mu ^{}\}2𝕀,`$ (16)
$`𝒜_{\text{conv}}`$ $`=`$ $`{\displaystyle \frac{k1}{\alpha \delta l}}{\displaystyle \underset{\mu }{}}{\displaystyle \frac{1}{1\alpha }}_\mu 𝕀,`$ (17)
$`𝒜_{\text{chem}}`$ $`=`$ $`{\displaystyle \frac{k}{\alpha }}{\displaystyle \underset{\mu }{}}{\displaystyle \frac{1}{1\alpha }}𝒞_\mu 𝕀.`$ (18)
Here the operators $`𝒜_\mu ^\pm ,_\mu ,𝒞_\mu `$ are defined according to the general definition (8) through their maps $`a_\mu ^\pm ,b_\mu ,c_\mu `$. In order to prove that the expectation value of the stochastic process whose generator is given by (13) really solves the Black-Scholes equation one uses Eq. (10) and the projection operator $`F_\lambda (\theta )=\theta _\lambda `$. One then immediately obtains

$`{\displaystyle \frac{}{t}}\theta _\lambda ={\displaystyle \frac{\theta _{\lambda +1}+\theta _{\lambda 1}2\theta _\lambda }{\delta l^2}}+(k1){\displaystyle \frac{\theta _{\lambda +1}\theta _\lambda }{\delta l}}k\theta _\lambda ,`$ (19)
which in the continuum limit $`\delta l0`$ converges towards the dimensionless Black-Scholes equation (6).
With these definitions the transition probability $`w(\theta ,\stackrel{~}{\theta })`$ becomes

$`w(\theta ,\stackrel{~}{\theta })`$ $`=`$ $`{\displaystyle \frac{1}{\alpha \delta l^2}}{\displaystyle \underset{\mu }{}}\delta (\theta ,a_\mu ^+(\stackrel{~}{\theta }))+\delta (\theta ,a_\mu ^{}(\stackrel{~}{\theta }))+{\displaystyle \frac{k1}{\alpha \delta l}}{\displaystyle \underset{\mu }{}}\delta (\theta ,b(\stackrel{~}{\theta }))`$ (20)
$`+{\displaystyle \frac{k}{\alpha }}{\displaystyle \underset{\mu }{}}\delta (\theta ,c(\stackrel{~}{\theta })).`$
Summarising the Black-Scholes Eq. (6) can be numerically solved by simulating the stochastic process

$`{\displaystyle \frac{}{t}}P(\theta ,t)`$ $`=`$ $`𝒜_{\text{BS}}P(\theta )={\displaystyle D\stackrel{~}{\theta }\left\{w(\theta ,\stackrel{~}{\theta })P(\stackrel{~}{\theta },t)w(\stackrel{~}{\theta },\theta )P(\theta ,t)\right\}}.`$ (21)
This is usually done with the help of the following algorithm:

* initialise $`\theta _\lambda `$ at $`t=0`$
* while $`t<t_{\text{end}}`$

+ compute random time step $`\tau `$ until the next jump occurs from an exponential distribution with mean value $`\tau =1/w_{\text{tot}}`$
    + apply one of the transitions $`a_\mu ^\pm ,b_\mu ,c_\mu `$ selected according to their relative probability

By repeating the above algorithm different realizations of $`\theta `$ are generated and thus $`\theta _\lambda `$ can be estimated from a sample of realizations.
This stochastic process has some interesting features. The parameter $`\alpha `$ can be used to control the size of the fluctuations. In earlier applications of the general theory presented above to balance equations of non-equilibrium statistical mechanics, the parameter $`\alpha `$ could be interpreted as the temperature of the physical system . The larger the parameter $`\alpha `$ the larger are the fluctuations, hence it is intuitively clear that, e.g. in a thermodynamic setting $`\alpha `$ plays essentially the role of the temperature.
From a numerical point of view the total transition probability $`w_{\text{tot}}`$ is constant. Thus the error made by taking a constant time step $`\tau =1/w_{\text{tot}}`$ in the numerical simulation vanishes as $`O(\delta t)`$. This significantly reduces the number of random variables needed. In addition the transition probabilities do not depend on the current state of the system. This makes the random selection of a transition very efficient, otherwise algorithms as discussed in  have to be used.
Fig. 1 shows a typical result of a simulation averaged over $`10`$ realizations of the stochastic process with $`\alpha =0.005`$. Naturally the simulation results fluctuate around the analytic solution, but these fluctuations are very small if one takes into account that only $`10`$ realizations are averaged. The reason for this of course is the parameter $`\alpha `$. The smaller $`\alpha `$, the smaller are the fluctuations of the process and the smaller is the number of realizations needed.
One drawback of this stochastic process is, that it is restricted to a uniform discretisation of the dimensionless variable $`\stackrel{~}{S}`$ from Eq. (6). As can be seen in Fig. 1 after transforming back to the original variables $`V`$ and $`S`$ this leads to an exponential distribution of the discrete values of the underlying, see Eqs. (5). Hence most of the points are at the left border of the integration region. Since the total transition rate performs as $`1/\delta l^3`$ the length of the time step during the simulation and thus the overall computing time depends strongly on the discretisation of the underlying. A numerically more efficient algorithm thus has to use a uniform discretization of the underlying $`S`$ which hence requires much less points. This results in a non-uniform discretization of the dimensionless variable $`\stackrel{~}{S}`$ in Eq. (6).

### 3.3 Stochastic process with non-uniform discretization

To enable a non-uniform discretization of the Black-Scholes equation (6) one uses an arbitrary distance $`\mathrm{\Delta }_\lambda `$ of two neighbouring points. The time evolution operator $`𝒜_{\text{chem}}`$ from the previous section does not depend on the discretization and is thus not altered. But for the diffusive and convective part of the Black-Scholes equation new stochastic processes $`𝒜_{\text{diff}}^{\text{n.u.}}`$ and $`𝒜_{\text{conv}}^{\text{n.u.}}`$ taking into account the non-uniform (n.u.) discretisation have to be defined. With the definitions

$`r_\mu ^+`$ $`=`$ $`\{\begin{array}{cc}\theta _\mu \hfill & \theta _\mu \alpha _1\frac{\theta _\mu }{\mathrm{\Delta }_{\mu 1}\mathrm{\Delta }_\mu }\hfill \\ \theta _{\mu +1}\hfill & \theta _{\mu +1}+\alpha _1\frac{2\theta _\mu \mathrm{\Delta }_{\mu +1}}{(\mathrm{\Delta }_{\mu +1}+\mathrm{\Delta }_\mu )\mathrm{\Delta }_\mu \mathrm{\Delta }_{\mu +1}}\hfill \end{array},`$ (22)
$`r_\mu ^{}`$ $`=`$ $`\{\begin{array}{cc}\theta _\mu \hfill & \theta _\mu \alpha _1\frac{\theta _\mu }{\mathrm{\Delta }_{\mu 1}\mathrm{\Delta }_\mu }\hfill \\ \theta _{\mu 1}\hfill & \theta _{\mu 1}+\alpha _1\frac{2\theta _\mu \mathrm{\Delta }_{\mu 2}}{(\mathrm{\Delta }_{\mu 1}+\mathrm{\Delta }_{\mu 2})\mathrm{\Delta }_{\mu 2}\mathrm{\Delta }_{\mu 1}}\hfill \end{array},`$ (23)
for the diffusive and

$`s_\mu `$ $`=`$ $`\{\begin{array}{cc}\theta _\mu \hfill & \theta _\mu \alpha _2\frac{\theta _\mu }{\mathrm{\Delta }_\mu }\hfill \\ \theta _{\mu 1}\hfill & \theta _{\mu 1}+\alpha _2\frac{\theta _\mu }{\mathrm{\Delta }_{\mu 1}}\hfill \end{array},`$ (24)
for the convective part one obtains the two new time evolution operators

$`𝒜_{\text{diff}}^{\text{n.u.}}`$ $`=`$ $`{\displaystyle \frac{1}{\alpha _1}}{\displaystyle \underset{\mu }{}}{\displaystyle \frac{1}{\text{Det}(r_\mu ^+)}}_\mu ^++{\displaystyle \frac{1}{\text{Det}(r_\mu ^{})}}_\mu ^{}2𝕀,`$ (25)
$`𝒜_{\text{conv}}^{\text{n.u.}}`$ $`=`$ $`{\displaystyle \frac{k1}{\alpha _2}}{\displaystyle \underset{\mu }{}}{\displaystyle \frac{1}{\text{Det}(s_\mu )}}𝒮_\mu 𝕀.`$ (26)
The total time evolution operator $`𝒜_{\text{BS}}^{\text{n.u.}}`$ is now given by

$$𝒜_{\text{BS}}^{\text{n.u.}}=𝒜_{\text{diff}}^{\text{n.u.}}+𝒜_{\text{conv}}^{\text{n.u.}}+𝒜_{\text{chem}}.$$
(27)
In order to prove that the expectation value $`\theta _\lambda `$ of the stochastic process above really solves the Black-Scholes equation one again introduces the projection operator $`F_\lambda `$ and makes use of the general theorem (10) to obtain

$`{\displaystyle \frac{}{t}}\theta _\lambda `$ $`=`$ $`{\displaystyle \frac{\frac{2\mathrm{\Delta }_{\lambda 1}}{\mathrm{\Delta }_\lambda +\mathrm{\Delta }_{\lambda 1}}\theta _{\lambda +1}+\frac{2\mathrm{\Delta }_\lambda }{\mathrm{\Delta }_\lambda +\mathrm{\Delta }_{\lambda 1}}\theta _{\lambda 1}2\theta _\lambda }{\mathrm{\Delta }_{\lambda 1}\mathrm{\Delta }_\lambda }}`$ (28)
$`+`$ $`(k1){\displaystyle \frac{\theta _{\lambda +1}\theta _\lambda }{\mathrm{\Delta }_{\lambda +1}}}k\theta _\lambda .`$
In the continuum limit the expectation value $`\theta _\lambda `$ of this stochastic process thus again solves the Black-Scholes equation, but now on a non-uniform grid. The transition probability $`w(\theta ,\stackrel{~}{\theta })`$ becomes

$`w(\theta ,\stackrel{~}{\theta })`$ $`=`$ $`{\displaystyle \frac{1}{\alpha _1}}{\displaystyle \underset{\mu }{}}\delta (\theta ,r_\mu ^+(\stackrel{~}{\theta }))+\delta (\theta ,r_\mu ^{}(\stackrel{~}{\theta }))+{\displaystyle \frac{k1}{\alpha _2}}{\displaystyle \underset{\mu }{}}\delta (\theta ,s(\stackrel{~}{\theta }))`$ (29)
$`+{\displaystyle \frac{k}{\alpha }}{\displaystyle \underset{\mu }{}}\delta (\theta ,c(\stackrel{~}{\theta })).`$
The parameters $`\alpha _i`$ are chosen such that the relative size of the transitions is smaller than $`1`$, hence they fulfil $`\alpha <1`$, $`\alpha _2/\mathrm{\Delta }_\mu <1`$ and $`2\alpha _1/\mathrm{\Delta }_\mu ^2<1`$ for every possible transition $`\mu `$. The smallest $`\alpha _i`$, which is $`\alpha _1`$, determines the scaling of the total transition probability

$$w=\frac{2n}{\alpha _1}+(k1)\frac{n}{\alpha _2}+\frac{k}{\alpha }.$$
(30)
Since $`\alpha _1O(\mathrm{\Delta }_\mu ^2)`$ and the total transition probability scales as $`O(n/\alpha _1)`$ the numerical effort to simulate the stochastic process scales according to $`O(n/\mathrm{\Delta }_\mu ^2)`$. This is the same dependence on the discretization as in the previous section with a stochastic process on a uniform grid. But since now the grid can be chosen appropriate one needs much less grid points and the algorithm is much faster.

### 3.4 Fast stochastic process with non-uniform discretization

The use of a non-uniform grid reduces the number of grid points to be used in a numerical simulation. But there is still a $`O(n/\mathrm{\Delta }_\mu ^2)O(1/\mathrm{\Delta }_\mu ^3)`$ dependence of the total transition rate which makes the algorithm slow. To get rid of this dependence the stochastic process corresponding to the diffusive part is again improved. In analogy to equilibrium Monte Carlo simulation of lattice systems , where one generates a new configuration of the whole lattice in a single sweep, one defines

$$𝒜_{\text{diff}}^{\text{sweep}}=\frac{1}{\alpha _1}\left[\frac{1}{\text{Det}(t^+)}𝒯^++\frac{1}{\text{Det}(t^{})}𝒯^{}2𝕀\right],$$
(31)
with the map $`t^\pm `$ given by

$$t^\pm =\underset{\mu }{}r_\mu ^\pm .$$
(32)
This process thus updates every $`\theta _\lambda `$ at once. The total time evolution operator is now given by

$$𝒜_{\text{BS}}^{\text{sweep}}=𝒜_{\text{diff}}^{\text{sweep}}+𝒜_{\text{conv}}^{\text{n.u.}}+𝒜_{\text{chem}}.$$
(33)
For the transition probability one obtains

$`w(\theta ,\stackrel{~}{\theta })={\displaystyle \frac{1}{\alpha _1}}\delta (\theta ,t^+(\stackrel{~}{\theta }))+{\displaystyle \frac{1}{\alpha _1}}\delta (\theta ,t^{}(\stackrel{~}{\theta }))`$ (34)
$`+`$ $`{\displaystyle \frac{k1}{\alpha _2}}{\displaystyle \underset{\mu }{}}\delta (\theta ,s(\stackrel{~}{\theta }))+{\displaystyle \frac{k}{\alpha }}{\displaystyle \underset{\mu }{}}\delta (\theta ,c(\stackrel{~}{\theta })).`$
Since $`\alpha _1`$ again scales as $`O(1/\mathrm{\Delta }_\mu ^2)`$ the total transition probability now also scales according to $`O(1/\mathrm{\Delta }_\mu ^2)`$. This of course makes the algorithm significantly faster.

### 3.5 Analysis of the algorithm

Fig. 2 shows three exemplary solutions of the stochastic process defined in (33) which corresponds to a non-uniform discretisation of the dimensionless underlying $`\stackrel{~}{S}`$. In addition to the expected random fluctuations around the exact solution, single realizations also show systematic deviations from the exact solution for a wide range of values of the underlying asset. These systematic and random fluctuations can be traced back to the two different types of stochastic processes entering Eq. (33). The generators $`𝒜_{\text{conv}}^{\text{n.u.}}`$ and $`𝒜_{\text{chem}}`$ describe single jump processes, where the jump process changes only one option value in asset time space, this of course causes the random fluctuations. In contrast the systematic deviations stem from $`𝒜_{\text{diff}}^{\text{sweep}}`$ which describes a jump process changing all option values at once.
The numerical estimation of the option price from a sample of realizations of the stochastic process is of course not affected by this behaviour of single realizations. By averaging about several realizations the estimated option price converges against the analytic solution, this is shown in Fig. 3.
Concerning the numerical root mean square error error

$$ϵ=\sqrt{\underset{k}{}\left[V(S_k,t)\widehat{V}(S_k,t)\right]^2},$$
(35)
of the Monte Carlo simulation of the Black–Scholes equations one has to distinguish two types of parameters entering the numerical algorithm derived in the previous section: the number of samples $`N`$ used to estimate the solution $`\widehat{V}`$ and the parameters $`\alpha _i`$ which control the size of the fluctuations. Figure 4 shows the root mean squared error $`ϵ`$ of a Monte Carlo simulation for three different values of the ratio $`\alpha _i/N`$. Different simulations of the same ratio were obtained e.g. by decreasing $`N`$ and all three parameters $`\alpha _i`$ by the same factor. This figure clearly shows that the error only depends on this ratio $`\alpha _i/N`$. In addition it can be seen that the error scales with the square root of this ratio. This of course is a dependence which is expected for a Monte Carlo simulation.

In order to assess the numerical performance of the proposed Monte Carlo method it is compared with the standard Monte Carlo approach based on the simulation of the dynamics of the underlying asset. The numerical error in a Monte Carlo solution has two sources. The systematic part of the error stems from the finite discretization of the stochastic differential equation in the case of the simulation of the geometric Brownian motion and from the finite discretization of the state space when one solves the master equation. In addition there is a random error from the averaging over different solutions which decays in both cases as $`N^{0.5}`$ (if $`\alpha `$ is fixed in the master equation formulation). Figure 5 shows the time required to achieve a given precision for the two methods in a parameter range where the systematic part of the numerical errors are comparable. As can be easily seen the solution of the master equation is significantly faster than the standard approach.

### 3.6 American options

Up to now only European options have been considered. How can the concept of mesoscopic Monte Carlo simulations be generalised for American options? A straightforward generalisation of the algorithm simulating the time evolution of the underlying asset is quite involved, since the early exercise condition has to be fulfilled. Whenever the option price falls below the payoff function the option is exercised and the option value is thus given by the payoff. Hence in order to assure that the early exercise condition is not violated sometime in the future a straightforward generalisation of this Monte Carlo simulation would have to keep track of all these points in asset time space which would make the algorithm very ineffective. Of course, there are more advanced techniques to generalise the Monte Carlo approach to American options, see  for an overview.
The proposed mesoscopic approach based on a master equation whose expectation value solves the original Black–Scholes equation can be generalised to the valuation of American options. Since the time direction in the dimensionless Eq. (6) has been reversed and the payoff function is used as initial condition the generalisation to American options is straightforward. One just has to simulate one time step according to the master equation (33) and average about the number of samples used to estimate the option price. Wherever this option price is below the payoff function it is replaced with the payoff and the next time step is simulated. Hence one assures that the early exercise condition is fulfilled.
For this approach to work it is critical that first the samples are averaged and that the early exercise condition is applied thereafter. Applying the early exercise condition before averaging would introduce a bias towards higher option prices. Every time the stochastic realization is by chance below the payoff function the option value is increased. But since the option value is never decreased this way it is intuitively clear that one would obtain higher option prices. Fig. 6 shows the result of such a simulation for an American call with exercise price 60, interest rate 0.07, volatility 0.2 and a continuous dividend yield of 0.10.

## 4 Conclusions

In real markets the elegance of the perfect hedge of the Black–Scholes approach to option pricing is generally lost. Already a very simple model of the dynamics of the underlying stock in terms of piecewise deterministic processes shows that the option price itself may be regarded as a stochastic variable. However the dynamics of the expectation value of this stochastic option price is governed by the Black–Scholes equation. This consideration is the motivation for the master equation approach to option pricing suggested in this paper.
The essential point of this approach is that the Black–Scholes equation is interpreted as the macroscopic equation of an underlying mesoscopic stochastic process for the stochastic option price variable. By using PDPs one can then easily proceed to a master equation formulation of the option pricing problem. In contrast to this the usually used microscopic approach describes the dynamics of the underlying asset by a stochastic differential equation.
The master equation formulation of the option pricing theory offers several advantages over the standard approach. This formulation provides a general setting in which also other kinds of jump processes may be easily embedded, without altering the character of the equations. One may also include additional stochastic processes for the volatility and the interest rate and then arrive at stochastic volatility and interest rate models . This of course would lead to different hedging strategies. One possibility to extend the proposed approach beyond the Black–Scholes equation, which is currently under investigation, is to use a hedging strategy proposed in
In addition the master equation formulation also allows for the use of advanced simulation algorithms. It was shown that it is possible to construct a master equation whose transition probability is constant in time and does not depend on the current state of the system. As shown in section 3.1 this allows for fast numerical algorithms for the computation of the option price.
Using the standard Black–Scholes equation as an example we have demonstrated how to construct such numerically efficient stochastic processes underlying the partial differential equation. As can be seen from the derivation of the general theory in section 3.1, it is clear that this approach is not restricted to the standard Black–Scholes equation but can be applied to generalisations of the latter.
The master equation formulation of the Black–Scholes equation is numerically about a factor of two faster than the standard Monte Carlo approach. This is mainly due to the fact that the proposed mesoscopic approach does not simulate sample trajectories of the underlying asset which, after averaging, results in the option price for the initial asset value at a given time. In contrast to this the mesoscopic approach works in the whole asset time space and generates sample option prices in the whole state space during one realization and is hence much faster.
Another advantage of the master equation formulation of the Black–Scholes equation is that it allows for a straightforward generalisation to price American options in the same framework by just comparing the average option price with the payoff function.
Summarising, the proposed master equation approach is located in between the standard Monte Carlo approach for the simulation of the underlying asset and the finite difference solution of the partial differential equation. This approach tries to combine the numerical efficiency of the solution of the partial differential equation with the advantages of the Monte Carlo approach. The main advantage of the Monte Carlo approach is that it can easily be applied to price options using also advanced stochastic models for the underlying.
Concluding, the mesoscopic approach we have formulated has been exploited to implement a master equation approach to option pricing. Of course, it will be of great interest to investigate in further work if features of real markets correspond to the additional freedom to choose the size of the fluctuations contained in the proposed master equation formulation. This will eventually lead to a better understanding of option pricing in real markets.


text:
# Interactions and correlations of particulate inclusions in a columnar phase

## I Introduction

Macroscopic particles suspended in an isotropic fluid - colloidal dispersions - have been the subject of intensive study for over a century colloid . Recently, a novel class of materials have been realised in which colloidal particles are dispersed in anisotropic fluids, typically thermotropic and lyotropic liquid crystals poulin ; zapot ; geetha . The combination of anisotropy and long-range order inherent in these mesophases leads to totally new kinds of interactions between the dispersed particles, which have no analogues in conventional colloidal dispersions brochard ; sriram ; terentjev ; poulinstark ; casimir .

The most important of these new interactions results from the elastic deformation induced by a particle in the liquid crystalline host. The overlap of the the elastic deformations of two particles produces an effective long-ranged interaction between them. This mechanism, which relies only on the mean distortions of the order parameter, implies that other sources of deformation like topological defects also interact with particles through long-range forces. Fluctuations of the order parameter together with the boundary conditions on the particle surface lead to another novel kind of Casimir-like interaction kardar , which can often be important near special points in the phase diagram. In particle-doped nematics, for instance, director fluctuations lead to a Casimir force which is stronger than the van der Waals force near the isotropic-to-nematic transition casimir . Interparticle interactions are further complicated by the inevitable presence of topological defects surrounding the particles poulinstark . These defects compensate the mismatch in the boundary conditions imposed on the order parameter at the particle surface and at infinity. To clarify with an example, a particle in a nematic with director-normal boundary conditions on its surface is equivalent to a radial hedgehog with unit positive topological “charge”, and is incompatible with the requirement of an uniform director configuration at infinity. The mismatch is relieved by the spontaneous appearance of a defect of unit negative “charge”, which is tightly bound to the particle. Such companion defects usually lead to rapidly decaying, repulsive, power law forces between the particles poulinstark . The interactions mentioned above, when combined with the multitude of interactions already present in conventional colloidal dispersions colloid , can drive the particles into spatially organised structures - chains poulinstark ; cladis , clusters raghu , and even crystalline arrangements nazaranko  have been observed for particles in nematic hosts.

As a result of these novel interactions, the static and dynamic properties of (with an obvious extension of terminology) liquid crystalline colloids are quite unusual. The change in the mechanical properties is particularly dramatic - particles added in volume fractions as low as several percent can totally alter the rheological behaviour of the host phase zapot ; geetha ; meeker . When the particles have a magnetic moment interesting ferro-liquid crystalline phases are formed ramosmag1 ; ramosmag2 .

Numerous theoretical studies have been done on these and related effects in nematic colloids nemcolloidref . Comparatively less work has been done on particle-doped lamellar phases pershan ; turner1 ; turner2 ; sens . It is therefore surprising that in spite of recent experiments ramosmag1 ; ramosmag2 , there has been no theoretical study of the properties of particle-doped columnar phases. In the following sections of this paper, therefore, as a first step towards a more detailed theory, we calculate the elastic field mediated interparticle interactions (section II), as well as the interactions of particles with point topological defects (section III). We conclude the paper with a critical discussion of our results and possible extensions of the present work.

## II Particles in a Columnar Hexagonal Phase

Let us recall that the columnar phase is macroscopically equivalent to a regular array of tubes, aligned in, say, the $`z`$ direction (directions in the $`x`$-$`y`$ plane are denoted as $``$) and free to slide along their axes. A cross-sectional slice normal to the tubes shows a two-dimensional crystal structure, which, for the case we study, is hexagonal. Mechanically, the phase is liquid-like in the $`z`$ direction, and solid-like and isotropic in the $``$ directions. Elastic deformations are described by a pair of “broken-symmetry” variables $`𝐮=(u_x,u_y)`$ corresponding to the local displacements of the tubes from their equilibrium position. To lowest order in gradients, the harmonic free energy density includes contributions from in-plane bulk and shear deformations, and out-of-plane column curvatures degennes :

$`_{hex}`$ $`=`$ $`{\displaystyle \frac{B}{2}}(_{}𝐮)^2+{\displaystyle \frac{K_3}{2}}\left(_z^2𝐮\right)^2`$ (1)
$`+`$ $`{\displaystyle \frac{C}{2}}[(_xu_x_yu_y)^2+(_yu_x+_xu_y)^2].`$
Here $`B`$ and $`C`$ are first-order elastic constants (with dimensions of energy density) corresponding to bulk and shear deformations respectively, while $`K_3`$ is a second-order elastic constant (with the dimension of force) corresponding to column curvature. An elastic description of this sort is valid at length scales large compared to the column spacing $`d`$. The colloidal particles are modelled as monodisperse spheres of radius $`a`$ and are described by a concentration field $`c(𝐫)=_\nu \delta (𝐫𝐫_\nu )`$, where $`𝐫_\nu `$ is the center-of-mass coordinate of the $`\nu `$th particle. A description of this sort is valid at length scales large compared to the particle radius $`a`$, when in effect, the particles can be treated as point masses. In this dilute limit, where the typical distance between particles is much larger than their size, all short-range interactions can be ignored, and the particles need to be assigned only an ideal gas entropy. As in the theory of point defects in three dimensional crystals, it is assumed that the particles behave as force dipoles (“centers of dilation or compression”) eshelby  and produce a body force proportional to $`c(𝐫)`$ at the point $`(𝐱_{},z)=𝐫`$. This is equivalent to a coupling of the form $`c(𝐫)E(𝐮)`$ in the free energy, where $`E(𝐮)=_{}𝐮`$ is the local dilation. It is not difficult to see that to linear order in $`𝐮`$ and lowest order in gradients, this is the only possible symmetry-allowed coupling of the concentration to the column displacement note1 . With these considerations, then, the free energy density (in units of $`k_BT`$) is

$$(c,𝐮)=_{hex}+c(𝐫)\left[ln\frac{c(𝐫)}{c_0}1\right]+\alpha c(𝐫)E(𝐮).$$
(2)
The second term in (2) accounts for the entropy of the particles whose mean concentration is $`c_0`$. $`\alpha `$ is a coupling constant,which on dimensional grounds must be of order $`Ba^3`$. Note that in our units $`B`$ has the dimensions of inverse volume and that $`\alpha `$ is dimensionless. The microscopic details of the boundary conditions on the particle surface are incorporated phenomenologically in the sign and magnitude of $`\alpha `$. For positive $`\alpha `$ the columns are repelled away from the particle surface leading to a local dilation; this situation is most appropriate for inert colloidal particles, whose size $`a`$ is appreciably larger than the column spacing $`d`$. For negative $`\alpha `$ the columns are are attracted towards the particle leading to a local compression; this situation may be appropriate for small biological inclusions, which are known to exert local inward forces on their surroundings turner1 ; safran . Henceforth, we shall assume $`\alpha `$ to be positive. It is worth emphasising that our theory is valid for separations large compared to both $`d`$ and $`a`$, which represent, respectively, the scales beyond which the elastic and concentration field descriptions are meaningful.

### II.1 Effective Interactions

Integrating out the displacement field $`𝐮`$ from $`F[c,𝐮]=d^3𝐫(c,𝐮)`$ we obtain an effective free energy

$$F_c=ln\left(𝒟𝐮e^{F[c,𝐮]}\right)$$
(3)
for the concentration fluctuations. Since our interest lies only in the mean displacements of the columns, the Gaussian integral in (3) can instead be replaced by the minimum of $`F[c,𝐮]`$ over $`𝐮`$. Setting the variation $`\delta F/\delta 𝐮`$ to zero, we obtain a pair of coupled linear differential equations for the components $`(u_x,u_y)`$. Transforming to Fourier space (ie multiplying by $`e^{i𝐪𝐫}`$ and integrating over $`d^3𝐫`$), and solving the resulting algebraic equation for the Fourier components, we have

$$𝐮_𝐪=\frac{i\alpha 𝐪_{}}{(B+C)q_{}^2+K_3q_z^4}c_𝐪.$$
(4)
This equation gives the (mean) displacement field for an imposed concentration fluctuation $`c_𝐪`$. Using it to eliminate the displacement field from (2), we obtain the effective free energy as

$$F_c=d^3𝐫c(𝐫)\left[ln\frac{c(𝐫)}{c_0}1\right]+\frac{1}{2}_𝐪c_𝐪V(𝐪)c_𝐪$$
(5)
where

$$V(𝐪)=\frac{\alpha ^2q_{}^2}{(B+C)q_{}^2+K_3q_z^4}$$
(6)
and we have used the shorthand $`_𝐪\frac{d^3𝐪}{(2\pi )^3}`$. The second term in (5) represents the effective interaction mediated by the column displacements. The kernel $`V(𝐪)`$ can be identified with the Fourier transform of the interparticle potential $`V(𝐫)`$, with the understanding that the correct physics is captured only for lengths large compared to the particle size. With that in mind, the interparticle potential is given by

$`V(𝐫)`$ $`=`$ $`{\displaystyle \frac{\alpha ^2\lambda _3^2}{K_3}}_{}^2G(𝐱_{},z)`$ (7)
$`=`$ $`{\displaystyle \frac{\alpha ^2\lambda _3^2}{K_3}}_{}^2{\displaystyle _0^{\mathrm{}}}𝑑t{\displaystyle \frac{e^{z^2/4\lambda _3x_{}\mathrm{cosh}t}}{(\lambda _3x_{}\mathrm{cosh}t)^{1/2}}}`$
where $`\lambda _3^2=K_3/(B+C)`$. The Fourier transform $`G(𝐱_{},z)`$ of the kernel $`(q_{}^2+\lambda _3^2q_z^4)^1`$ can be obtained explicitly in terms of modified Bessel functions prost , but we prefer the integral representation given above.

The potential $`V(𝐫)`$ presents several interesting features. It is non-central (not a function of $`r`$ only) and radially non-monotonic. This last property prevents us from inferring the sign of the radial component of the force from that of the potential. Therefore, the attractive or repulsive character of the interaction cannot be correlated with the sign of the potential, and the forces have to be considered explicitly. However, in the $`z`$ and $``$ directions the potential is monotonic and decays with a power law. Near the $`z`$-axis (for $`z^24\lambda _3x_{}`$) the potential vanishes as $`|z|^5`$, leading to a rapidly damped and repulsive force along the columns.Near the $`x`$-$`y`$ plane (for $`z^24\lambda _3x_{}`$) the potential decays more slowly as $`x_{}^{5/2}`$, implying a repulsive force in the $``$ direction. In a general direction the nature of the interaction is best visualised through the plot of the equipotentials of $`V(𝐫)`$ shown in Fig.1. From the plot it is clear that the radial component of the force is positive near the $`z`$-axis and the $`x`$-$`y`$ plane, and negative in the region bounded by the two paraboloids $`z^2\lambda x_{}`$. Broadly speaking, the force felt by a particle at $`(𝐱_{},z)`$ due to another particle at the origin varies from repulsive, to attractive, to repulsive again as the polar angle $`\theta =\mathrm{tan}^1(z/x_{})`$ varies from $`0`$ to $`\pi /2`$. Since the potential is even in $`z`$, the sequence repeats as the angle decreases from $`\pi `$ to $`\pi /2`$. A qualitatively similar angular dependence is also found in the quadrupolar forces between particles in a nematic sriram ; the distance dependence, though, is totally different.

Given the nature of the interparticle forces described above, we might expect an attraction induced aggregation of the particles. However, repulsive forces due to the companion topological defects mentioned in the introduction may dominate at short distances, and arrest the aggregation. In that case, it is likely that particles will spontaneously organise in chains oblique to the columns. This process will be dominated by Brownian motion, and not the elastic interaction, which though long-ranged, decays too rapidly to affect the kinetics of aggregation.

The correction to the interparticle potential due to a finite density of particles may be calculated quite easily within a mean field approximation cl . Augmenting the free energy $`F_c`$ with an external potential $`V_{ext}(𝐫)`$ which couples to the concentration, the equilibrium concentration is obtained by demanding $`\delta F_c/\delta c=0`$. This gives

$$c(𝐫)=c_0\mathrm{exp}\left[V_{ext}(𝐫)d^3𝐫^{}c(𝐫^{})V(𝐫𝐫^{})\right].$$
(8)
The mean (or effective) potential is defined by $`c(𝐫)=c_0e^{U(𝐫)}`$. Using this to eliminate the density from the previous equation and identifying the external potential with that of a fixed particle at the origin, we obtain the following self-consistent equation for $`U(𝐫)`$:

$$U(𝐫)=V(𝐫)+c_0d^3𝐫^{}V(𝐫𝐫^{})\mathrm{exp}\left[U(𝐫^{})\right].$$
(9)
Assuming that $`U1`$ we can linearise the above integral equation. Neglecting the term corresponding to the background concentration and Fourier transforming we find

$$U(𝐪)=\frac{V(𝐪)}{1+c_0V(𝐪)}=\frac{\alpha ^2q_{}^2}{(B+Cc_0\alpha ^2)q_{}^2+K_3q_z^4}.$$
(10)
The mean potential is of the same functional form as the ”bare” interparticle potential, but with a “renormalised” characteristic length $`\lambda _3^2(\alpha )=K_3/(B+Cc_0\alpha ^2)`$. Recalling that a similar calculation with a Coulomb potential $`V(𝐪)q^2`$ leads to charge screening cl , we must conclude that, despite the long-ranged interaction, many-particle effects are quite modest in our case.

### II.2 Correlation Functions

In the absence of interactions the particles in our model have only configurational entropy and their two-point correlations are like those of an ideal gas: $`<\delta c(𝐫)\delta c(𝐫^{})>\delta (𝐫𝐫^{})`$. However, elastic interactions lead to quite complex correlations, as we show below. Defining the dimensionless structure factor $`S(𝐪)`$ as the Fourier transform of $`\frac{1}{c_0}<\delta c(𝐫)\delta c(𝐫^{})>`$ we find from the effective free energy in (5) that

$`S(𝐪)`$ $`=`$ $`{\displaystyle \frac{1}{1+c_0V(𝐪)}}`$ (11)
$`=`$ $`1+{\displaystyle \frac{c_0\alpha ^2q_{}^2}{(B+Cc_0\alpha ^2)q_{}^2+K_3q_z^4}}.`$
The Fourier transform of the pair correlation function $`h(𝐪)=S(𝐪)1`$ is proportional to $`U(𝐪)`$, and therefore the pair correlations themselves are well described by the contour plot of Fig.1 with all the signs reversed. The correlations are negative along the $`z`$ and $``$ axes, indicating a reduction in the particle concentration compared to the mean value of $`c_0`$, while in the region bounded by the paraboloids $`z^2\lambda _3x_{}`$ the correlations are positive, indicating an enhancement of the concentration. This is broadly consistent with our earlier picture of particle chains oblique to the columns. It may be possible to measure the particle structure factor by small-angle X-ray scattering, and the form of $`S(𝐪)`$ above would then provide a test of the particle-column coupling assumed in our theory.

While the column displacements induce particle correlations, the particle aggregates in turn modify the column fluctuations. In particle-doped lamellar phases, this feedback mechanism has been shown to change the usual Caillé caille  form of the structure factor sens . To calculate the column displacement correlations in the presence of particles, we first define the displacement field such that $`_iu_j=0`$ for a mean concentration $`c_0`$. This requires that we modify the particle-column coupling to $`\delta c(𝐫)E(𝐮)`$. Next, we define an effective free energy $`F_𝐮`$ for the column fluctuations (in a way similar to the definition of $`F_c`$) by integrating out the concentration field. Carrying out the integration by steepest descent we obtain the effective free energy

$$F_𝐮=d^3𝐫\left[_{hex}c_0\left(\alpha _{}𝐮e^{\alpha _{}𝐮}\right)\right].$$
(12)
This shows that the coefficients of all symmetry-allowed terms which contain $`_{}𝐮`$ are “renormalised” by the presence of particles. Within our harmonic description the only relevant term, obtained by expanding the exponential term to quadratic order, is $`\alpha ^2(_{}𝐮)^2`$. Retaining only this term in (12) and defining the longitudinal and transverse parts of $`𝐮_𝐪`$ as $`q^2𝐮_L=𝐪(𝐪𝐮_𝐪)`$ and $`𝐮_T=𝐮_𝐪𝐮_L`$ we obtain the $`u_{L(T)}u_{L(T)}`$ correlations as

$$|𝐮_T|^2=\frac{1}{Cq_{}^2+K_3q_z^4},$$
(13)
$$|𝐮_L|^2=\frac{1}{(B+Cc_0\alpha ^2)q_{}^2+K_3q_z^4}.$$
(14)
The transverse correlations are identical to those of the pure columnar phase, which indicates that there are no corrections to $`C`$ and $`K_3`$ due to the particles. The longitudinal correlations while identical in form, have a modified bulk modulus $`B_c=Bc_0\alpha ^2`$. Particles thus produce a softening of the columnar phase.

## III Particle-defect interactions

### III.1 Point defects in the columnar phase

Topological defects are common in liquid crystalline phases, being easily produced under conditions of processing such as shear flows.Much of the interesting macroscopic behaviour of liquid crystalline colloids comes from the interaction of the colloidal particles with the topological defects of the host zapot ; geetha ; meeker . For example, the solid-like response of particle-doped cholesteric zapot  and smectic geetha  phases is due to the formation of a stable network of defects. The colloidal particles play the crucial role of stabilising this network, which otherwise would anneal away. It is likely that elastic interactions do not play a dominant role in such processes, and that short-range forces, not captured in a coarse-grained description such as ours, are more important. Nevertheless, elastic interactions between particles and topological defects are known to have important consequences kosevich , and so in what follows, we calculate this interaction for a point defect configuration in the columnar phase. The statistical mechanics of thermally activated point defects jain ; nelson  may be modified in the presence of particles and a knowledge of particle-defect interactions could be potentially useful in a study of this phenomenon. An interesting duality (see below) between the columnar and lamellar phases has been noted by de Gennes and Prost degennes  and explored in detail by Selinger and Bruinsma bruinsma . This suggests that it may be of interest to compare particle-defect interactions in the lamellar and columnar phases and we conclude this section with such a comparison.

The defect configuration we consider is a single column ending prost ; bruinsma , illustrated in Fig. 2. This is not unduly restrictive, since any linear defect which consists of column endings (eg. a transverse edge dislocation) can be constructed out of a collection of such point defects prost . The ease with which a column can break depends on what constitutes the “column”, and therefore, on what the particular realisation of the columnar phase is. In polymeric columnar phases, the columns consist of long chains arranged end-to-end along their length. Here, thermally activated column scission can occur only at those places where one end of a chain meets that of another. (Polymeric liquid crystals are found predominantly in the nematic and isotropic phases; however, longer and stiffer biopolymers like DNA can also positionally order into columnar phases polymerhex ). Since columns in discotic and micellar phases are composed of weakly bonded molecules, thermally activated scission can take place anywhere along the length of the columns. In all of these examples, the column scission energy is finite and thus a state of thermal equilibrium at non-zero temperature will always contains a finite number of such defects. The ends of a scissioned column may either separate along the $`z`$-axis to form a string of vacancies, or may translate past each other in the $``$ direction to form a string of interstitials jain . Such strings will begin from and end on the column extremities, the lower of which, by convention, is called a “tail”(Fig. 2) and the upper a “head” jain . At high temperatures the entropy of wandering of the string can overcome the energy due to its line tension, leading to a string proliferation, and producing isolated “heads” or “tails” prost ; frey . Signatures of such defects have been observed in recent experiments endpoint . An isolated semi-infinite string running along the $`z`$-axis with its terminus at the origin imposes the constraint that bruinsma

$$_{}𝐮=m\frac{\delta (𝐱_{})}{\rho _\mathrm{𝟎}}\theta (z),$$
(15)
where the topological “charge” $`m`$ is $`+1`$ ($`1`$) for a tail (head), and $`\rho _0`$ is the areal density of the columns. Clearly this is nothing but a line of centers of dilation and is easily treated in the framework of the previous section. Formally, a column terminus at the origin is identical to a line of point particles distributed with the concentration

$$c(𝐫)=m\frac{(B+C)}{\alpha \rho _0}\delta (𝐱_{})\theta (z).$$
(16)
Combining this with the expression in (4) for the displacement produced by a distribution of particles, the displacement due to a column ending is obtained as note2

$$𝐮_𝐪=\frac{i𝐪_{}}{q_{}^2+\lambda _3^2q_z^4}\frac{m}{\rho _0iq_z}.$$
(17)
The displacement for a collection of column endings distributed with the topological “charge” density $`m(𝐫)`$ follows from the superposition of the individual displacements. Stated mathematically, this means that the total displacement is

$$𝐮_𝐪=\frac{i𝐪_{}}{q_{}^2+\lambda _3^2q_z^4}\frac{m_𝐪}{\rho _0iq_z},$$
(18)
with $`m_𝐪`$ being the Fourier transform of $`m(𝐫)`$. To calculate the particle-defect interaction, the displacement field is decomposed as $`𝐮_𝐪=𝐮_𝐪(c_𝐪)+𝐮_𝐪(m_𝐪)`$, where the first part (due to the particles) is obtained from (4) and the second part (due to the defects) is obtained from (18). The free energy then separates out as $`F=F_c+F_{mm}+F_{cm}`$, with

$$F_{mm}=\frac{1}{2}_𝐪\frac{q_{}^2}{\rho _0^2q_z^2}\frac{m_𝐪m_𝐪}{q_{}^2+\lambda _3^2q_z^4}\frac{K_3}{\lambda _3^2}$$
(19)
and

$$F_{cm}=_𝐪\frac{i\alpha q_{}^2}{\rho _0q_z}\frac{m_𝐪c_𝐪}{q_{}^2+\lambda _3^2q_z^4}.$$
(20)
The Fourier transforms of the defect-defect and particle-defect interactions can now be read off as the kernels of $`m_𝐪m_𝐪`$ and $`c_𝐪m_𝐪`$ in (19) and (20) respectively. The transformation to real space can be simplified by noting that

$$_{}^2G(𝐱_{},z)=\lambda _3^2_z^4G(𝐱_{},z).$$
(21)
Using this property, $`q_{}^2`$ can be replaced by $`\lambda _3^2q_z^4`$ in the numerators of the kernels in (19) and (20). Alternatively, a conjugate Green’s function can be introduced as in bruinsma  and nelson . With either of the methods, the final form of the interaction kernels turn out to be

$$V_{mm}(𝐪)=\frac{K_3}{\rho _0^2}\frac{q_z^2}{q_{}^2+\lambda _3^2q_z^4},$$
(22)
and

$$V_{cm}(𝐪)=\frac{\alpha \lambda _3^2}{\rho _0}\frac{iq_z^3}{q_{}^2+\lambda _3^2q_z^4}.$$
(23)
Apart from the line tension terms which we ignore, the defect-defect interaction is identical to the result of bruinsma  and has been discussed in detail there note3 . The particle-defect interaction can be represented in real space as

$$V_{cm}(𝐫)=\frac{m\alpha \lambda _3^2}{\rho _0}_z^3G(𝐱_{},z).$$
(24)
As might have been anticipated, the particle-defect potential is non-central and radially non-monotonic, and is odd in the $`z`$-coordinate. The behaviour in the $`z`$ and $``$ directions is simple. Along the columns (for $`z^24\lambda _3x_{}`$) the potential decays as sgn($`z`$)$`|z|^4`$ indicating a repulsive force for $`z>0`$ and an attractive force for $`z<0`$. In the $`x`$-$`y`$ plane, the $``$ component of the force vanishes by symmetry and the $`z`$ component goes to zero as $`x_{}^{5/2}`$. The detailed form of the interaction for positive $`z`$ is shown in the equipotential plot of Fig.3. The general conclusion is that independent of the sign of the defect, the particles are attracted to the regions of compression, and repelled away from the regions of dilation that the defect alone would have produced, see Fig.2. Driven by the attractive component of the interaction, and aided by Brownian motion, a cloud of particles should collect around a column ending, a situation reminiscent of the Cottrell clouds of vacancies and interstitials around dislocations in metals kosevich . Normally, such clouds decrease the dislocation mobility kosevich  and the same can be assumed to happen here. Since defect motion is one of the main mechanisms of stress relaxation in ordered phases, less mobile defect strings may lead to an increase in the viscosity of the columnar phase.

It is important to note that expression (18) can be used to obtain the interaction of a particle with any linear defect which can be represented as a collection of column endings. For example, the interaction of a transverse edge dislocation with a particle can be obtained from (18) by setting $`m(𝐫)\delta (y)\delta (z)`$ and repeating the steps leading to (20).

### III.2 Line defects in the lamellar phase

There is an interesting line-plane duality degennes ; bruinsma  between defect-free columnar and lamellar phases which is apparent in a pictorial representation: a set of lines representing the columns (aligned, say, along the $`z`$-axis) when rotated by $`\pi /2`$ can equally well be viewed as a set of planes representing the lamellar layers (now stacked normal to the $`z`$-axis). Mechanically, the lamellar phase is a 2d liquid and a 1d solid, while the columnar phase is a 1d liquid and a 2d solid. Though this duality is lost in the presence of an arbitrary ensemble of defects (for example screw dislocations), nonetheless there is one important exception. It is easy to see from Fig. 2 that an endpoint defect (the terminus of a line) in the columnar phase corresponds to an edge dislocation (the terminus of a plane) in the lamellar phase. Going beyond pictures, there are close connections in the energetics of these “dual”defects, which have been explored in bruinsma . This prompts us to ask if similar connections may be found in the interactions of endpoints and edge dislocations with particles. To answer this question, we first examine the particle-dislocation interaction for an arbitrary distribution of dislocations, and then specialise to the case of an edge dislocation.

The free energy density of a particle-doped lamellar phase where the particles couple to the local layer dilation has been obtained in turner1 . Excluding the entropic terms, in our notation it is

$$F=\frac{\overline{B}}{2}d^3𝐫\left[(_zu)^2+\lambda ^2(_{}^2u)^2+\overline{\alpha }c(𝐫)_zu\right],$$
(25)
where $`u`$ is the local displacement of the layers from their equilibrium position, $`B`$ is the bulk modulus, $`K=B\lambda ^2`$ is a second order elastic constant corresponding to layer curvature, and the positive constant $`\overline{\alpha }`$ is the strength of the particle-layer coupling. In the presence of a dislocation density $`𝐛(𝐫)`$ the distortion $`𝐯=u`$ acquires a singular transverse part determined by pershan

$$\times 𝐯=𝐛(𝐫).$$
(26)
The total distortion in Fourier space is given by cl

$$𝐯_𝐪(𝐛_𝐪)=\frac{i𝐪\times 𝐛_𝐪}{q^2}+i𝐪h_𝐪,$$
(27)
where $`h_𝐪`$, the longitudinal part of the distortion, is obtained from the equation of equilibrium $`\delta F/\delta u=0`$. The solution

$$h_𝐪=\frac{\overline{\alpha }iq_zc_𝐪}{q_z^2+\lambda ^2q_{}^4}\frac{q_z(1\lambda ^2q_{}^2)}{q_z^2+\lambda ^2q_{}^4}\frac{\left(\widehat{𝐳}𝐪\times 𝐛_𝐪\right)}{q^2}$$
(28)
is a sum of the particle and defect contributions. When the distortion in (27) is inserted into the (Fourier transformed) free energy (28), the latter separates into three pieces: the first contains only the dislocation density cl , the second only the concentrations turner1 , and the third is a cross-term which represents the particle-dislocation interaction. Displaying only this last term, we have

$$F_{cb}=\overline{B}_𝐪\frac{\overline{\alpha }\lambda ^2q_{}^2}{q_z^2+\lambda ^2q_{}^4}(\widehat{𝐳}𝐪\times 𝐛_𝐪)c_𝐪.$$
(29)
This gives the interaction energy of an arbitrary distribution of particles with an arbitrary number of dislocations. Two important special cases are linear edge and screw dislocation of infinite extent. The dislocation density for a screw dislocation of infinite extent and Burgers vector $`b\widehat{𝐳}`$ is $`𝐛(𝐫)=b\delta (x)\delta (y)\widehat{𝐳}`$, which implies that $`𝐛_𝐪=2\pi b\delta (q_z)\widehat{𝐳}`$. Inserting this into (29) we find that there is no elastic interaction between a particle and the screw component of a dislocation within a linear elastic theory. For the case of interest, an infinite edge dislocation running along the $`y`$-axis with Burgers vector $`b\widehat{𝐳}`$, $`𝐛(𝐫)=b\delta (x)\delta (z)\widehat{𝐲}`$, and $`𝐛_𝐪=2\pi b\delta (q_z)\widehat{𝐳}`$, and so together with (29) the particle-edge dislocation interaction kernel is obtained as

$$V_{cb}(𝐪)=\frac{\overline{B}\overline{\alpha }\lambda ^2q_{}^2}{q_z^2+\lambda ^2q_{}^4}2\pi biq_x\delta (q_y).$$
(30)
Transforming to real space, the interaction potential is

$$V_{cb}(𝐫)=\frac{\overline{B}\overline{\alpha }b}{8}\frac{x}{\sqrt{\pi \lambda |z|^3}}\mathrm{exp}\left(\frac{x^2}{4\lambda |z|}\right).$$
(31)
The interaction is long-ranged along the $`z`$-axis (for $`x^24\lambda |z|`$) and exponentially damped along $`x`$-axis (for $`x^24\lambda |z|`$). By symmetry, the $`z`$-component of the force vanishes on the $`z`$-axis, while both the $`x`$ and $`z`$ components are exponentially small near the $`x`$-axis. Irrespective of the sign of $`b`$, the particles are attracted to the compression regions that the dislocation alone would have produced, and repelled away from the dilation regions, see Fig.(4). To compare this interaction with its dual in the columnar phase, we first note that the interaction is characterised by the parabolas $`x^2=4\lambda |z|`$ which occur frequently in problems related to dislocations in lamellar phases degennes . The particle-defect potential in the columnar phase is characterised instead by the paraboloids $`z^2\lambda _3x_{}`$. In both cases the maximum elastic deformation is confined to the regions within these surfaces, compare Fig.(2) and Fig.(4). The force is long-ranged along the solid-like directions (the $``$ direction in the columnar phase, and the $`z`$ direction in the lamellar phase) and short-ranged in the liquid-like directions. There is no interaction between a particle and an edge dislocation when both have the same $`z`$ coordinate, a feature also seen in the interaction between two edge dislocations degennes . In contrast, the particle-defect interaction in the columnar phase is strongly repulsive in the $`z`$ direction. The formation of Cottrell clouds, and the consequent reduction in the dislocation mobility should be expected in the lamellar phase as well.

## IV Conclusion

To summarise, we have calculated the elastic interactions of colloidal particles in a columnar hexagonal phase by modelling them as centers of dilation. The interaction is long-ranged and non-central, and has both attractive and repulsive parts. We have shown that this interaction leads to correlations between the particle positions, and to a reduction in the bulk modulus $`B`$ of the columnar phase. We have also shown how particles interact elastically with defects in columnar and lamellar phases. In both phases, particles aggregate around the defects driven by the attractive component of the interaction and reduce the defect mobility. One likely consequence of this is an increase in the viscosity of the phase.

We mention below several of the directions in which the present work may be extended. First, the elastic interactions of particles of more general shapes like cylinders and discs should be examined. Such inclusions can no longer be modelled as force dipoles, and force distributions containing higher multipoles will be needed eshelby . The interaction of plate-like particles would have a direct relevance to the physics of clay-block copolymer composites groenewold . Second, the effects of boundaries should be investigated. Free boundaries require the introduction of image sources, which will have long-range interactions with particles in the bulk. Thus, boundary effects can be expected to strongly influence all processes in the bulk, and in particular, the formation of particle aggregates. Third, the nature of the topological defects surrounding a particle should be examined. Our coarse-grained description cannot predict the defect configuration (though useful information can be gained if we assume it) and a more microscopic description would be better suited for the purpose.

As a final, somewhat speculative remark, we should indicate that the present calculation may be of some relevance in biological systems where inclusions may be present in tissues of columnar hexagonal structure bouligand . Biological inclusions generate local inward forces turner1 ; safran , and a minimal description of such inclusions could be in terms of the centers of compression (corresponding to a negative $`\alpha `$) in our theory. Similar elastic interactions between proteins in a lamellar phase turner2 , and cells in an isotropically elastic tissue have already been investigated safran .

###### Acknowledgements.

The author would like to thank Professor Sriram Ramaswamy for suggesting the problem, useful discussions, and a critical reading of the manuscript, and Professor Chandan Dasgupta for helpful suggestions. Partial support from Hindustan Lever Limited is gratefully acknowledged.