\section{}

\section{Introduction}

Rare-earth compounds provide an exceptional frame to study magnetic interactions. 4f and 5f shells are incomplete giving often magnetic moments. Magnetic ordering is controlled by the interplay between exchange interactions and crystal field coupling \cite{Canfield06}. Some rare earth based materials are superconducting and magnetic \cite{Canfield06,Crespo09}, with phase diagrams showing a mutual relationship between metamagnetic and superconducting transitions \cite{Galvis12}. These phase diagrams have brought considerable interest \cite{Canfield98, Buzdin85, Aoki01, Petit10}.\\

To further advance, it is useful to characterize magnetism in simple metallic systems with a crystal structure which is easy
to handle. Of interest, in particular for first principles calculations, are the monopnictide compounds. They crystallize in a NaCl-like structure, what makes them suitable to perform theoretical analysis. Therefore these compounds have attracted interest for testing magnetic interaction theories \cite{Mullen74, Pawlicki88} and more recently due to their semiconducting properties and practical applications \cite{Bhardwaj11}. A structural transition into a CsCl structure appears with pressure \cite{Bhardwaj11, Svane01}. The electronic phase diagram of an entire range of rare earth monopnictides \cite{Petit10} has been predicted using \textit{ab-initio} calculations. Experiments and theory have addressed several systems also in the particular case of Ho-monopnictides, such as HoP \cite{Fischer85_1}, HoAs \cite{Schmid85} or pure and Y doped HoSb \cite{Kim76, Hessel80, Jensen80, Taub74}.\\

HoBi is a metallic compound. Its lattice parameter is of $a=0.62291(3) nm$ \cite{Pearson} and its magnetic behavior stems from the interaction between the Ho ions. Ho$^{3+}$ has total angular momentum of J=8 and the magnetic behavior is expected to be correspondingly rich. Previous studies on HoBi single crystals, grown from the melting of polycrystalline samples, find at low temperatures \cite{Ott84} an antiferromagnetic transition at 5.7 K based on thermal expansion and x-ray measurements, and a structural transition into a pseudotetragonal cubic structure with an axial ratio $c/a$ greater than unity. Neutron diffraction studies \cite{Fischer85} performed in polycrystalline samples find a second order transition into a fcc type II antiferromagnet. Recent studies use density functional theory to predict structural an electronic properties of this compound \cite{Coban10}. Bhajanker {et al} using also theoretical calculations have predicted a structural transition under pressure at 26 GPa into a CsCl-like structure for HoBi \cite{Bhajanker12}.
None of the studies we found about HoBi addresses the evolution of its magnetic phases with both field and temperature.\\

In this paper we describe the synthesis and characterization of high quality single crystals of HoBi. We study in particular temperature and magnetic field dependence of the magnetization with the aim to build a H-T phase diagram at low temperatures. We find typical features of a second order antiferromagnetic transition and a number of field-induced transitions. Values found for the effective Ho-magnetic moment and saturation magnetization agree with previously reported ones  \cite{Ott84}.\\

\section{Material and Methods}

\begin{figure}[ht!]
	\begin{center}
		\includegraphics[width=0.4\textwidth,clip]{./Figure1.eps}
	\end{center}
	\caption{(a)HoBi crystal grown from excess of Bi flux. Each of the small squares in the background corresponds to one millimeter. (b) Crystallographic structure of HoBi where blue (big) and red (small) spheres represent Ho and Bi ions respectively.}
	\label{crystal}
\end{figure}

High quality single crystals of HoBi were grown from Bismuth flux \cite{Growth_1, Growth_2} using a 10\% atomic percentage of Ho and 90\% of Bi. Pure elements were placed into an alumina crucible and heated from room temperature to $1000\,^{\circ}\mathrm{C}$ in three hours, cooled in another three hours to $900\,^{\circ}\mathrm{C}$ and slowly cooled in 104 hours to $400\,^{\circ}\mathrm{C}$. As a result we obtained crystals like the one in figure \ref{crystal} with sizes between 0.2 mm and 1 mm. Crystals show a cubic-like shape according to their $NaCl$ structure and are easily cleaved along the main crystallographic axes.\\

Magnetization measurements were performed in a Quantum Design Magnetic Properties Measurement System (MPMS) in zero field cooling and with the external field applied parallel to the [100] direction. M(H) curves were measuring stabilizing the temperature and going up and down in field. Since practically no hysteresis was found we only display here measurements when increasing the magnetic field. Regarding M(T) curves we stabilized the field and we went up in temperature. For resistance measurements we use a Quantum Design Physical Properties Measurement System (PPMS) and the four probe method with platinum wire contacts glued with silver epoxy. HoBi is an air sensitive material and in order to reduce exposure we did not polish or shaped the sample. This condition makes it hard to give an accurate resistivity value due to the relevance of the geometric factor in this calculation. Zero field specific heat was measured between 2 K and 20 K using adiabatic microcalorimetry again with a PPMS system.\\

\section{Results and Discussion}

\begin{figure}[ht!]
	\begin{center}
		\includegraphics[width=0.9\columnwidth,clip]{./Figure2.eps}
	\end{center}
	\caption{In (a) we show the temperature dependence of the resistance normalized to its ambient temperature value at zero magnetic field. Current is applied along the a axis of the cubic structure. In (b) we show the temperature dependence of the inverse susceptibility with a magnetic field of 0.1 T applied along the a axis. Blue line is a fit to the expression given in the text between 55 K and 350 K. Inset: Zero field specific heat shows a jump at the magnetic transition.}
	\label{quality}
\end{figure}

The temperature dependence of the resistance of HoBi normalized to its ambient temperature value is shown in Fig.\ref{quality}.a from 300 K to 1.8 K. The resistance linearly drops down to 50 K, and then saturates. Below 5.7 K, at the magnetic transition temperature, we observe a significant drop. The RRR (Residual Resistance Ratio) is 126 between 300 K and 1.8 K.\\

Fig.\ref{quality}.b shows the inverse susceptibility at a magnetic field of 0.1 T. We find a linear Curie-Weiss-like behavior until HoBi enters the magnetically ordered phase at 5.7 K. This transition is consistent with the drop in resistance and therefore we can associate the latter to a significant decrease of magnetic scattering. Magnetic ordering increases the electron mean free path \cite{Canfield95}.\\

We can fit the susceptibility as a function of temperature using
	$$ \chi = \frac{C}{T-\Theta_C}, $$
and obtain between 55 K and 350 K a Curie-Weiss temperature of $\Theta_C$=-0.26 K. The effective Bohr magneton is extracted from
	$$C=\frac{N\mu_B^2}{3Vk_B}\ p^2 ,$$
giving a value of p=11.1. The free Ho$^{3+}$ ion has a Bohr magneton per Ho of 10.6. We find a value slightly bigger than the expected one but in the same order. Below the transition the sample presents the typical temperature dependence of the susceptibility for an antiferromagnet.\\

The specific heat (inset of Fig.\ref{quality}.b) shows a clear lambda like anomaly with a peak at the magnetic transition temperature again consistent with the transition we observed in susceptibility and resistance measurements and with a dependence very similar to the one found for HoSb in ref. \cite{Taub74}.\\

\begin{figure}[hb!]
	\begin{center}
		\includegraphics[width=0.9\columnwidth]{./Figure3.eps}
	\end{center}
	\caption{(a) M/H for different fields show antiferromagnetic behavior below the transition temperature for the lower fields. Between H=2 T and H=4 T different metamagnetic states are reached at the lowest temperatures and for field above 4.00 T no transition is seen. (b) We show $d(MT/H)/dT$ to better highlight the transition temperatures \cite{Fisher62}. Data for different temperatures have been shifted by 0.001 cm$^3$/mol. (c) Detail of $d(MT/H)/dT$ for some fields to show the types of features found in this derivatives and the criteria used for the selection of the transition temperatures. The H=2 T and the H=2.25 T data are shifted by 5 (cm$^3$/mol) and -5 (cm$^3$/mol) respectively for better visualization.}
	\label{trans_temp}
\end{figure}

From these plots we can identify the antiferromagnetic transition of HoBi at zero field at 5.7 K. When we apply a magnetic field along the [100] direction, we observe that magnetism in HoBi evolves strongly with the magnetic field. We find several field-induced metamagnetic states. Fig.\ref{trans_temp}.a shows the evolution of temperature dependent M/H for several fields between 0.10 T and 5.00 T. For increasing fields the transition is broadened and moved to lower temperatures. The new metamagnetic states emerge for fields higher than 1.50 T. For fields above 4 T, metamagnetic transitions are suppressed. Applying a magnetic field above 5.00 T leads to a saturated paramagnetic state where all the moments are aligned.\\

In order to extract the transition temperatures from the $M/H$ data the temperature derivative of $MT/H$ \cite{Fisher62} is plotted in Fig.\ref{trans_temp}.b. Maxima and shoulders of this derivative were taken as the transition temperatures \cite{Tomy95}. Low field plots show sharp maxima that can be associated to the entrance into the antiferromagnetic state. As the field is increased first a high temperature and then a low temperature shoulder (both associated with the metamagnetic states) emerge and the peak becomes smaller and broader. For higher fields the peak disappears completely and only the two shoulders remain. Fig.\ref{trans_temp}.c shows three representative curves of the series containing
 the different features we associate with the transition temperatures.\\

\begin{figure}[htb!]
	\begin{center}
		\includegraphics[width=0.9\columnwidth,clip]{./Figure4.eps}
	\end{center}
	\caption{(a) Magnetization measurements with the data shifted by 1$\mu_B$ between different temperatures. The low temperature data show a cascade of up to five transitions that disappear as the temperature is increased. (b) We plot the derivative of the magnetization with the magnetic field as a function of the magnetic field. Peaks correspond to transitions, providing the phase diagram discussed in the next figure. Derivative data have been shifted by $2 \mu_B/Ho T$ for the different temperatures.}
	\label{trans_field}
\end{figure}

We have also measured magnetization as a function of the magnetic field at different fixed temperatures up to 5.50 T. Previous measurements in Ref. \cite{Ott84} find significant hysteresis in most transitions. Here we find practically no hysteresis. We ascribe this difference to the crystal growth method, which, in our case, should lead to more homogeneous samples. In our measurements, at saturation, magnetization at low temperatures reaches the 10$\mu_B$ expected from Ho$^3+$ ions within the mass error. Fig.\ref{trans_field}.a show the magnetization for temperatures between 1.85 K and 6.00 K (above the zero field transition temperature). The low temperature data show a cascade of field-induced transitions that disappear as the temperature is increased until no transition is seen at 6.00 K.\\

Metamagnetic transitions in M(H) data are easily determined from the derivative dM/dH (Fig.\ref{trans_field}b). Five clean transitions can be seen as maxima in 1.85 K and 2.00 K data. For 2.50 K only three maxima can be seen but a broad feature emerges around 2 T and its trace can be followed up to at least 4.00 K. Between 3.00 K and 5.00 K two maxima can be followed, while only one remains at 5.50 K. At 6.00 K magnetization is a straight line and all transitions have disappeared.\\

\section{Conclusions}

\begin{figure}[htb!]
	\begin{center}
		\includegraphics[width=0.9\columnwidth]{./Figure5.eps}
	\end{center}
	\caption{(a) In the H-T phase diagram several magnetic states can be differentiated. For low fields HoBi is paramagnetic (PM) down to $T_N$ when it enters the antiferromagnetic (AF) phase. When higher external fields are applied up to five metamagnetic (MM) states appear. (b) Detail of the lower temperatures phase diagram. Lines are guides to the eye.}
	\label{phase_diagram}
\end{figure}

With the data from $dM/dH$ and $d(MT/H)/dT$ we obtain a magnetic phase diagram shown in Fig.\ref{phase_diagram}. For the region above 3.00 K three different areas are clearly differentiated: the paramagnetic (PM) state above $T_N$, the antiferromagnetic (AF) below this transition and a metamagnetic state that is labeled as MM1. Below 3.00 K we can differentiate four more metamagnetic states labeled MM2-5.\\

In Fig.\ref{phase_diagram}.b we show a zoom over the low temperature region of the phase diagram. We observe five phases, which emerge at low temperatures from two magnetic phases. 2.5 K seems to be a characteristic temperature, where the magnetic order becomes more complex and field dependent. Fields higher than 5 T force all magnetic moments to point in the [100] direction, leading the system into an induced-ferromagnetic state.\\

Reference \cite{Ott84} studied HoBi single crystals finding a cascade of six metamagnetic transitions at 1.5 K in a squeezed sample. These transitions appear for fields lower than the ones we find here. From Maxwell relation $\left(\left( \partial S/\partial P \right)_T = - \left(\partial V/\partial T \right)_P\right)$, and the positive thermal expansion found in ref. \cite{Ott84}, we deduce that pressure and stress should give lower transition temperatures.\\

To summarize, we have grown and characterized single crystals of HoBi using the excess flux method. Our crystals have a very high RRR, evidencing a low concentration of defects and high electron mean free path. We have built a phase diagram showing a number of stable metamagnetic phases and their evolution with increasing field and temperature. The magnetic phase diagram of HoBi is amazingly rich and deserves further study to understand the spin structure and series of transitions until saturation.\\

\section{Acknowledgments}
We thank discussions with J.V. Alvarez. This work was supported by the Spanish MINECO (Consolider Ingenio Molecular Nanoscience CSD2007-00010 program, FIS2011-23488 and MAT2011-27470-C02-02), by the Comunidad de Madrid through program Nanobiomagnet and by NanoSc-COST program. We also acknowledge Banco Santander. Work at the Ames Laboratory was supported by the US Department of Energy, Basic Energy Sciences, Division of Materials Sciences and Engineering under Contract No. DE-AC02-07CH11358.\\

\section{References}
\bibliographystyle{model1a-num-names}


\section{Introduction}

Recent initiatives for expanding human presence in the solar system, including further exploration of the Moon and Mars, have generated interest in understanding planetary crater formation from a retrograde rocket \cite{metzger1,metzger2,metzger4,metzger5,metzger3,Kuang}. This problem essentially amounts to a jet of gas impinging on a granular bed from above. In addition to its practical application, this fundamental two-phase flow (gas and granular material) comprises an interesting problem in it's own right. Depending on the strength and composition of the jet and the characteristics of the granular material, the crater formation can occur in different ways (e.g. surface erosion of particles, penetration of the gas into the granular material, bulk shear of the granular bed) \cite{metzger1}. However, we will restrict ourselves primarily to considering so-called viscous erosion \cite{Kok}, in which the jet lifts or rolls the top few layers of grains, and a scour hole forms beneath the jet. This is the easiest regime to reach in that it takes the smallest jet velocity to accomplish. Additionally, due to a pulse when the jet is switched on, the initial formation of the craters in our experiments contains a component of eruption from gas diffusing into the interstices. The strength of this effect depends strongly on the velocity of the jet. In our experiments, we used a high-speed camera to observe the crater formation in a quasi-2D environment,  characterizing the process by the time-evolution of the shape, size, and location of the crater, as we varied control parameters of the system. During this process, we observed an unexpected symmetry breaking: under certain conditions, the centralized crater becomes unstable, and the crater moves to the left or to the right. Figure~\ref{symbreakstrong} shows frames from videos, where symmetric craters break symmetry in two different modes which we denote ``strong'' and ``weak'', which will be the focus of the bulk of this paper. We will attempt to give a coherent picture of the cratering process as a whole, including the symmetry-breaking effect, which is consistent with the data and gives physical insights into the process, at least on a phenomenological level.

\begin{figure*}[h!tbp]
\centering
\includegraphics[trim=10mm 0mm 10mm 0mm, width=\textwidth]{symbreakframesweak_small.eps}
\includegraphics[trim=8mm 0mm 8mm 15mm, width=\textwidth]{symbreakframes_small.eps}
\caption{Still frames of a weak, leftward (TOP) and a strong, rightward (BOTTOM) symmetry breaking process (for both sets of frames: left to right, at 3, 5, 7, and 15 seconds). Weak symmetry breaking typically occurs with higher jet heights and lower pressures, while strong symmetry breaking typically occurs with lower jets and higher pressures. The weak case shown has a jet starting height of $h=75$~mm and pressure of $P=7$~psi (see text for details of experimental setup). Throughout the process, the crater retains its shape fairly well and debris is ejected from both sides of the crater. For the strong case shown, with $h=35$~mm and $P=10$~psi, the debris is ejected more vertically, leading to a strong coupling between the two sides of the crater. As the crater shifts to the right, debris is ejected only from the right side of the crater, and a circular flow pattern emerges (counterclockwise, in images shown here).}
\label{symbreakstrong}
\end{figure*}

\section{Experimental Setup}

The primary apparatus used in our experiments (Fig. \ref{SetupPic}) consists of a thin, transparent acrylic box
with inner dimensions approximately $25\times 50\times 1$ cm. The box was partially
filled with various granular materials, each of which possesses different characteristics which clearly
affected the dynamics of the cratering process, including an earth-source red sand, as well as lunar
and Martian soil simulants manufactured by Orbital Technologies (ORBITEC). A thin metal pipe (22 cm long, with inner diameter of
0.44 cm) was inserted through a hole in the top of the box. The pipe was connected to a pressure-controlled
tank of nitrogen gas through a flexible tube: 40 cm with inner diameter 0.89 cm, then a valve, then 37 cm with inner diameter 0.38 cm to the pressure sensor and nitrogen tank. The total length of the path between the pressure sensor and the nozzle is $\sim 1$~m. We varied $P$, the pressure driving the gas flow between 2 and 10 psi, $h$, the starting height of the jet above the granular material between 35~mm and 75~mm, and the type of material used, i.e., the three materials discussed above. Using high-speed videos at 100 frames per second, we characterized the cratering process by measuring $D$, the depth of the crater from the initial surface as a function of time, $w$, the width and general shape of the crater as it evolved, and $A_x$, the horizontal displacement of the bottom of the crater with respect to center, which served as a measurement of the horizontal asymmetry.

\begin{figure}[h!tbp]
\centering
\includegraphics[width=0.9\columnwidth]{SetupSketch.eps}
\caption{Sketch of the experimental apparatus}
\label{SetupPic}
\end{figure}

Note that this apparatus is different from the ones used in previous experiments by Metzger et al. \cite{metzger1,metzger2,metzger4,metzger5}, in that it is quasi-two-dimensional instead of a three-dimensional half-space (where the jet driving the crater formation is placed against a transparent side wall of a container, effectively taking a cross-section of a 3D cratering experiment). These two approaches are both chosen for their obvious advantages over the fully 3D case in terms of taking measurements. By making either simplification, one obviously changes the system. For the primary crater formation, the quasi-2D system exhibits behavior which is very similar to the 3D half-space.

The material properties were as follows: the sand consisted of roughly monodisperse (about 1 mm diameter) but rough particles with moderate macroscopic friction (internal angle of friction, $\delta=37$ degrees), the lunar simulant, JSC-1, was very frictional ($\delta=45$ degrees) and had a broad distribution of particle sizes (0.01 - 2 mm) (see \cite{mckay,metzger5}). The
Mars simulant, JSC Mars-1, was less frictional then the lunar simulant ($\delta=36$ degrees) and also had a broad size distribution
(0.02 - 1 mm) \cite{allen}.

\section{Experimental Results}

\subsection{Symmetric Crater Evolution}

A general qualitative description of the process by which cratering occurred is as follows: the impinging jet of gas created an elevated pressure and an upward flow at the bottom of the crater. The flow entrained a relatively thin layer of granular material along the sides and bottom of the crater, and drove this material out of the crater. As one might expect, the smaller, lighter particles tended to become entrained in the gas flow much more easily than heavier particles. Hence, in materials with a large distribution of sizes, such as the lunar or Martian simulants, the flow of debris exiting the crater typically contained a higher concentration of smaller particles. We note that the regime of gas flow in our experiments was one where the gas had enough kinetic energy to move the particles, but not enough kinetic energy to penetrate the granular material in a way that caused eruption.

We consider first what is probably the most obvious measurement in a cratering process, namely, the depth of the crater, $D(t)$. In all the data,
the depth of the crater was found to grow logarithmically with time during the initial stages:

\begin{equation}
D(t) =a \ln{(b t +1)}
\label{eq:depthlog}
\end{equation}

This basic form governs the depth evolution for a substantial period of time after the jet is turned on, and it is applicable across a wide range of parameter values for jet speed, jet starting height, and material used, as long as the digging process is one of viscous erosion. The logarithmic growth in these quasi-2D experiments is similar to the observations by Metzger et al. \cite{metzger1} in a 3D half-space. Logarithmic growth is observed in cratering studies with widely varying experimental conditions \cite{metzger1,metzger2,rajaratnam,mazurek,varas,kellay}, which suggests an underlying physical principle which manifests itself irrespective of the details.

The formation of the crater was not simply determined by ejection of material, but rather, involved competition between digging and refilling processes. As the crater grew deeper and wider, ejected debris had a greater chance of being recirculated back into the crater, by either falling directly in the crater or by avalanching in from an outer crater which formed in some cases and is discussed below. After a period of logarithmic growth, the crater eventually approached an asymptotic depth, $D_{\infty}$,
where the ejection of matter by driving gas is balanced by a reflux of grains back into the crater. In these experiments, this saturation point seemed to depend on the box geometry, including its size. We return later to this issue, where we suggest an evolution equation which captures the initial logarithmic growth and the asymptotic depth. Characteristic data is shown in Fig. \ref{saturation}.

\begin{figure}[h!tb]
\centering
\includegraphics[width=\columnwidth]{symcratergrowth.eps}
\caption{Above is data demonstrating the asymptotic crater depth. The crater grows logarithmically in time
for the first 10 seconds (top). However, observing for a longer time, we clearly see the depth approach an asymptote, $D_\infty$ (bottom). This data is from a cratering experiment with red sand in which there was no breaking of the horizontal symmetry, as will be discussed later. We later discuss the dashed red line below in the context of a growth model which allows for saturation.}
\label{saturation}
\end{figure}

Eq.~\ref{eq:depthlog} requires a characteristic length scale ($a$), and a characteristic time scale ($1/b$). An important issue concerns how the dependence of $a$ and $b$ depend on other system parameters. In a 3D half-space experiment \cite{metzger2}, Metzger, et al, investigated the dependence of these parameters on the starting height of the nozzle
($h$), the kinetic energy density or momentum flux of the gas ($\rho v^2$), and the inner diameter of the nozzle ($d$). These authors found that $a$ scales linearly with $h$ and is unaffected by $\rho v^2$ or $d$, while $b$ scales linearly with $\rho v^2$, $h^{-3}$, and $d$. Our measurements agree with this: if we assume that the jet velocity scales linearly with the driving gas pressure, $P$, implying that the flow in the connecting line is close to laminar, our results are consistent with the Metzger et al. measurements.

\begin{figure}[h!tbp]
\centering
\includegraphics[width=0.85\columnwidth]{lunarcrater.eps}
\includegraphics[trim=0mm 0mm 0mm 5mm, width=0.85\columnwidth]{marscrater.eps}
\caption{Images taken during crater formation using lunar simulant (TOP) and Mars simulant (BOTTOM). For the lunar simulant, the crater shape is approximately a hyperbolic cosine (dashed line). The shape remains very consistent as the crater grows: depth and width both grow logarithmically in time with extremely similar time scales. For the Mars simulant, as well as for the sand used for the bulk of this paper (shown in Fig.~\ref{symbreakstrong}), reduced friction gives rise to inner and outer crater.}
\label{lunarcrater}
\end{figure}

The basic shape of the crater was most clearly affected by the frictional properties of the granular material, as in \cite{metzger5}. Less frictional materials, such as the sand and Mars simulant, Figs.~\ref{lunarcrater} (bottom) and Fig.~\ref{symbreakstrong}, formed craters with two clear sections: an inner crater, which is directly below the impinging jet and had rounded shape (discussed below), and the outer crater, which connected the outer edge of the inner crater along planes to the surface at approximately the angle of repose of the granular material. Highly frictional materials, i.e. the lunar simulant, did not exhibit a significant outer crater, at least not for time scales of these experiments. It seems likely that the outer crater is formed due to Coulomb failure when the slope of the side wall becomes steep. The more frictional materials can support a nearly vertical side wall with some added help from the pressure of the jet (see Figure \ref{lunarcrater}), whereas the less frictional materials cannot.

Also of interest is the shape of the inner crater, which we characterize by the aspect ratio $\Gamma$, defined as the width-to-depth ratio of the inner crater.
As shown in Fig.~\ref{aspectratio}, the starting height of the jet has a strong influence on the fundamental geometry of the crater. In fact, Fig.~\ref{aspectratio} shows that $\Gamma$ scales linearly with the jet starting height, with at best a weak dependence on air speed/pressure. This dependence of $\Gamma$ on jet height can be seen in the top-left and bottom-left images of Fig.~\ref{symbreakstrong}: the top case has $h=75$~mm and $\Gamma\approx 2$, and the bottom case has $h=35$~mm and $\Gamma\approx 1$, consistent with Fig.~\ref{aspectratio}. We can approximate $\Gamma$ with the following form:
\begin{equation}
\Gamma = (h - h_0)/H_0,
\label{geom}
\end{equation}
where $h_0$ is about twice the tube inner diameter, and $H_0 \simeq 2.2$.

\begin{figure}[ht]
\centering
\includegraphics[scale=0.9]{aspectratio.eps}
\caption{Measurements of the aspect ratio for craters in red sand after 1 second: the jet height sets the crater's shape. The aspect ratio, $\Gamma$ which is measured as the width-to-depth ratio of the inner crater, is linearly dependent on jet height and only very weakly dependent on driving pressure.}
\label{aspectratio}
\end{figure}

\subsection{Horizontal Symmetry Breaking}

As previously stated, for certain conditions, the cratering process became horizontally asymmetric, with the crater moving to one side or the other. We observe two qualitatively different modes of horizontal symmetry breaking in the cratering process, shown in Fig.~\ref{symbreakstrong}: the top shows weak symmetry breaking, which is characterized by broad craters and slow transitions, and the bottom shows strong symmetry breaking, which is characterized by narrower craters and rapid transitions. This phenomenon was observed in the red sand and in the Martian simulant, but not in the Lunar simulant, even for very strong jets. This difference was presumably due to the frictional properties of the materials (that is, symmetry breaking requires a less frictional material). The following results concerning symmetry breaking are all from experiments using the red sand.

To quantify the symmetry breaking, we measure the depth, $D$, and horizontal deviation from the center, $A_x$, of the lowest point of the crater. Figure~\ref{weakdata} shows these quantities for a weak symmetry-breaking case. The logarithmic growth of depth is not substantially disrupted by the symmetry breaking, and at the same time the horizontal position of the crater changes gradually, i.e. quasi-linearly.

\begin{figure}[h!]
\centering
\includegraphics[trim=0mm 0mm 0mm 0mm, width=0.9\columnwidth] {symbreakweakdepth.eps}
\includegraphics[trim=1mm 0mm 0mm 6mm, width=0.9\columnwidth] {symbreakweakhorz.eps}
\caption{Symmetry breaking, with $P=6$~psi and $h=35$~mm, for weak symmetry breaking, with depth (top) and asymmetry (bottom) plotted: points are experimental data, dotted lines show logarithmic and linear fits respectively, and the horizontal asymmetry, $A_x$, is the deviation from the center position. The lowest point of the crater moves quasi-linearly as it drifts from symmetric to asymmetric. At the end of the linear transition, the horizontal position asymptotes to its fully asymmetric state, where one edge of the inner crater is aligned with the jet.}
\label{weakdata}
\end{figure}

In contrast, for the strong symmetry-breaking case, the logarithmic growth of the depth is severely disrupted when the symmetry is broken, and the horizontal position changes much more rapidly. Fig.~\ref{strongdata} shows the depth and horizontal positions of craters for two separate experiments for the same initial conditions. Both data sets start along the same path, where $D(t)$ follows a logarithmic function and $A_x$ is zero. The red data follows this trajectory much longer than the blue data. Each set spontaneously breaks from the logarithmic function at a different time. At this time, the depth growth is at least temporarily interrupted and the asymmetry, $A_x$, begins to grow, following a logarithmic evolution curve, relatively smoothly at first, and then with fluctuations imposed. During the latter stages of evolution, the depths for the two different cases approach each other, and likewise the horizontal positions of the two data sets approach each other. In fact, the rather complex structures for depth and asymmetry can be reproducible if the symmetry breaking occurs at the same time after initiating the flow. For instance, Fig.~\ref{strongdata2} shows two experiments with the same initial conditions which happen to break symmetry at the same time. Note the consistency in the behavior, both in the symmetric growth period (logarithmic depth growth) and in the symmetry breaking period (quasi-logarithmic horizontal change).

\begin{figure}[th]
\centering
\includegraphics[trim=0mm 0mm 0mm 0mm, width=0.9\columnwidth] {symbreakdiffdepth.eps}
\includegraphics[trim=0mm 0mm 0mm 5mm, width=0.9\columnwidth] {symbreakdiffhorz.eps}
\caption{Symmetry breaking data for strong symmetry breaking, depth (top) and asymmetry (bottom): this figure shows two experiments (plotted in red and blue) with identical initial conditions ($P=6$~psi, $h=35$~mm). The symmetry breaking occurs at different times for the two different experiments. The symmetry breaking shows no preferential behavior as far as when it will occur or which
side it will move toward.}
\label{strongdata}
\end{figure}

\begin{figure}[th]
\centering
\includegraphics[width=0.9\columnwidth] {symbreaksamedepth.eps}
\includegraphics[trim=0mm 0mm 0mm 6mm, width=0.9\columnwidth] {symbreaksamehorz.eps}
\caption{Data for strong symmetry breaking, depth (top) and asymmetry (bottom): this figure shows two experiments (plotted in red and green) with identical initial conditions ($P=7$~psi, $h=35$~mm) such that symmetry breaking occurs at the same time after initiation of the air flow. Note the similarity between the two experiments for both the depth, $D$ and asymmetry, $A_x$.}
\label{strongdata2}
\end{figure}

There is an additional insight that can be made by observing different types of digging and the corresponding crater geometry. As discussed previously, a larger initial height of the nozzle yields a broader crater, and a lower nozzle yields a narrower crater. Most importantly, as mentioned previously, material leaves more vertically from narrower craters than for broad craters. The direction in which granular material is ejected has the effect of essentially coupling or decoupling the debris leaving the two sides of the crater and then returning. For broad craters, it is extremely unlikely that debris leaving one side of the crater ends up on the other side. However, for very narrow craters, a significant portion of the debris exiting one side can be directed to the opposite side, which can strongly affect horizontal asymmetries. That is, if there is a narrow crater where the coupling between the two sides is extremely strong, once the digging becomes slightly more prominent on one side, a substantial portion of that debris is dumped back into the other side. This effect appears to drive the strong symmetry-breaking case.

We identify strong symmetry-breaking as the case where the sidewalls of the inner crater range from nearly vertical to overhanging and a circular-like flow pattern emerges: the flow tends to go down one edge of the crater and upward and out the other, throwing debris back into the downward flow and into the other side of the crater. Likewise, we define weak symmetry-breaking as the case where the sidewalls are not as steep and there is not a circular-like pattern, but where there is a clear preferential direction to the digging, albeit not as severe. With these definitions in mind, Fig. \ref{symbreakdata} indicates when each type of symmetry-breaking occurs as a function of the starting height of the nozzle and the driving pressure. In summary, symmetry breaking does not occur for very large $h$ and/or small $P$. When it occurs, smaller $h$ is necessary for a given $P$ if strong symmetry breaking is to occur.

\begin{figure}[h!tbp]
\centering
\includegraphics[width=0.9\columnwidth]{symbreakdata.eps}
\caption{Plot shows where different kinds of symmetry breaking tend to occur as a function of the starting nozzle height ($h$) and the driving pressure for the gas ($P$). Dashed lines are meant to suggest different regions where each type of cratering occurs, but further investigation would be required to adequately map out these regions.}
\label{symbreakdata}
\end{figure}

This symmetry breaking is strongly suggestive of a pitchfork bifurcation associated with the the asymmetry, $A_x$, in the horizontal direction. In such a bifurcation, an increase in the relevant control parameter causes a single stable fixed point to become unstable to two stable symmetric fixed points. Such a bifurcation could be either forward or backward, where the latter is characterized by hysteresis and a discontinuous jump in the steady state value of $A_x$ as the control parameter is raised above the transition.

First we look at the continuity of the of $A_x$ during symmetry breaking. In the strong symmetry breaking case, it appears discontinuous (see Figs. \ref{strongdata} and \ref{strongdata2}), while in the weak symmetry breaking case, it seems much slower and less abrupt (see Fig. \ref{weakdata}). This suggests that perhaps the strong symmetry breaking is backward and the weak symmetry breaking is forward.

We can also look for hysteresis to confirm these classifications. By operating the system in an approximately adiabatic mode, where the pressure is varied slowly in small increments so the depth stays approximately at its asymptotic value, one can clearly observe hysteresis at a jet height (and crater shape) characteristic of strong symmetry breaking processes. Figure~\ref{hysteresis} shows data from an example of this procedure, with a jet height of 35 mm. As the pressure is slowly increased, the crater grows deeper. At some point, the system breaks symmetry. Once the symmetry breaking occurs, the pressure is held constant until the system stabilizes. Then, the pressure is slowly decreased to a level well below the value for the original symmetry breaking. However, the amount of asymmetry remains effectively constant.

For weak symmetry breaking (i.e. with a jet height of 75 mm and a wide crater shape), no sudden symmetry breaking occurs in an adiabatic mode of operation. Instead, the crater might shift a very small amount, but it is always possible to reduce the pressure again and have the system remain stable and relatively symmetric. Note that even a symmetric cratering process is innately hysteretic and irreversible (i.e., if a significant crater is created, the system can never return to a crater-less state simply by changing the gas flow in some way, such as varying the pressure or jet height). In this way, hysteresis exists in both types of symmetry breaking processes.

\begin{figure}[h!tbp]
\centering
\includegraphics[width=\columnwidth]{hysteresis.eps}
\caption{Characteristic data for an ``adiabatic'' process (see text) with a jet height $h=35$~mm, where the pressure is varied sufficiently slowly between 3 and 6~psi, such that the crater evolution follows a succession of nearly steady states. The data set begins just before the onset of symmetry breaking. At 6~psi, the crater breaks symmetry and runs to one side, as the crater deepens slightly. Video shows this to be a strong symmetry breaking transition, and the system parameters ($P=6$~psi, $h=35$~mm) are consistent with the strong regime in Fig.~\ref{symbreakdata}. After the crater has stabilized, the pressure is slowly decreased to a value which is below the initial onset of symmetry breaking. The crater remains fully asymmetric during this process, showing no signs of returning to the middle.}
\label{hysteresis}
\end{figure}

\section{Discussion and Conclusions}

\subsection{Analysis of Logarithmic Growth}

Given the ubiquity of a logarithmic form for cratering from fluid flows,
it is useful to consider mathematical models, i.e. one or more differential equations, which would give rise to such a form. For instance cratering experiments including those of Metzger et al., the present experiments, and even experiments with a disturbance coming from below a submerged granular bed \cite{varas}, modeling a volcanic eruption, exhibit logarithmic growth.
At the simplest level, one might write the following evolution equation for $D$:
\begin{equation}
\frac{dD}{dt}=a b \exp{(-D/a)}.
\label{depthexp}
\end{equation}
This equation contains both a length scale, $a$, and a time scale, $b$. The length scale may well be tied to the spreading of the gas plume as it approaches the granular bed (see discussion above). The gas speed (determined by the pressure driving the flow) and $a$ are then sufficient to determine a time scale. In fact, a somewhat different form for the evolution
of $D$ is also useful. Defining $u = \dot{D}$, yields an evolution equation
\begin{equation}
\frac{du}{dt} = -u^2/a
\label{ueq}
\end{equation}
which does not contain an explicit time scale. Integrating to obtain $u(t)$, yields
\begin{equation}
u = \frac {dD}{dt} = u_o a/(a + u_o t)
\end{equation}
where $u_o$ is the initial value of $dD/dt$, which has units of velocity. Identifying $u_o$ with the gas jet velocity seems reasonable, in the absence of additional information. The logarithmic form for $D(t)$ then follows as
\begin{equation}
D(t)= a \ln[1 + (u_o t/a)].
\label{newlogform}
\end{equation}
We note that the cutoff in the logarithmic growth can easily be modeled by:
\begin{equation}
\frac{dD}{dt}=a b \left[\exp{(-D/a)}-\exp{(-D_\infty /a)}\right],
\label{modifieddiffeq}
\end{equation}
where $D_\infty$ is the saturation depth. This equation leads to the following time dependence:
\begin{equation}
D(t) =a \ln{\left(e^{\frac{D_\infty }{a}}-(e^\frac{D_\infty }{a}-1)e^{-bt \exp {(-D_\infty /a)}}\right)}.
\label{modifieddepth}
\end{equation}
For short times, we recover the previous form:
\begin{equation}
D(t) \approx a \ln{(bt+1)}.
\end{equation}
For long times, the system approaches an asymptotic value of
\begin{equation}
D(t\rightarrow\infty) =D_\infty .
\end{equation}

Revisiting the data in Fig.~\ref{saturation}, the functional form of Eq.~\ref{modifieddepth} is shown to be an acceptable fit to the crater depth in the logarithmic and saturation regimes. These results are indicated by the thin red line.

\subsection{Conclusions}

The present experiments have explored the nature of cratering caused by impinging gas jets on quasi-2D samples of various granular materials, including lunar and Martian simulants and terrestrial sand. We observe initially logarithmic grow of the crater depth. We also observe what we believe is a previously unreported symmetry breaking instability. Assuming that a similar effect occurs in 3D, this symmetry breaking could be of extreme practical concern for the design of a spacecraft which needs to land on or take off from a granular surface. A lander design based on symmetric crater evolution could lead to unexpected catastrophic results if a symmetry breaking instability occurs and the crater moves laterally, disrupts the surface where the lander should make contact with the ground, or if the exhaust from the rocket, along with entrained material, is redirected towards the lander. For instance, an imbalance or redirection of the thrust or exhaust stream might exert enough torque to flip it over. We also note that previous work \cite{metzger6,metzger7} has suggested that planetary landings, especially on less frictional granular surfaces (i.e., Mars), may have other potential failure modes.

An interesting analogue may occur in experiments involving water droplets hitting a bed of sand covered by a thin layer of water. In the experiments of Kellay \cite{kellay} droplets created a circularly symmetric disturbance for a 3D sand bed, and at a certain drop frequency, the circular symmetry was broken, leading to a logarithmically evolving asymmetric pattern.

\begin{acknowledgements}

This work was supported by ORBITEC (contract \# OTC-GS-02381), subcontracted from USAF (contract \# NNX09 CF72P)
and by NASA contract. We particularly thank Dr. Phil Metzger for input on various aspects of the problem of crater formation.

\end{acknowledgements}


\section{Introduction}

Recently, the both schools of Takemura and Takayama have developed a quite interesting minimization method called holonomic gradient descent method(HGD). It utilizes Gr\"{o}bner basis in the ring of differential operator with rational coefficients. Gr\"{o}bner basis in the differential operators  plays a central role in deriving some differential equations called a Pfaffian system for optimization. HGD works by a mixed use of Pfaffian system and an iterative optimization method. It has been successfully applied to several maximum likelihood estimation (MLE) problems, which have been intractable in the past. For example, HGD  solve numerically the MLE problems for the von Mises-Fisher distribution and the Fisher-Bingham distribution on the sphere (see,  Sei et al.(2013)  and Nakayama et al.(2011)). Furthermore, the method has also been applied to the evaluation of the exact distribution function of the largest root of a Wishart matrix, and it is still rapidly expanding the area of applications(see, Hashiguchi et al.(2013)).
On the other hand, in statistical models, it is not rare that parameters are constrained and therefore the MLE problem with constraints has been surely one of fundamental topics in statistics. In this paper, we develop HGD for MLE problems with constraints, which we call the constrained holonomic gradient descent(CHGD).
The key of CHGD is to separate the process into (A) updating of new parameter values by Newton-Raphson method with penalty function and (B) solving a Pfaffian system.

\section{Constrained Optimization Problem}

We consider the following the constrained optimization problem.
\begin{eqnarray}
(P)&min&\ f(x) \nonumber \\
&s.t.& g_i(x) \leq 0,\ h_i(x)=0
\end{eqnarray}
where $i=1,...,m$, $j=1,...,l$
and $f,g_i,h_j : R^n \rightarrow R$ are all assumed to be continuously differentiable function. $g_i(x)$ is an equality constraint function and $h_j(x)$ is an inequality constraint function. In this paper, the objective function $f(x)$ is assumed to be holonomic.
We call the interior region defined by the constraint functions {\it the feasible region}.

\subsection{Penalty Function Method}

A penalty function method replaces a constrained optimization problem by a series of unconstrained problems.
It is performed by adding a term to the objective function that consists of a penalty parameter $\rho$ and a measure of violation of the constraints. In our simulation, we use {\it the exact penalty function method}. The definition of the exact penalty function is given by (see Yabe (2006)).
\begin{eqnarray}
P(x; \rho):= f(x) + \rho\{\sum_{i=1}^m |g_i(x)| + \sum_{j=1}^l max(0,h_j(x))\}, \ \rho>0
\end{eqnarray}

\section{Holonomic descent method}

Assume that we seek the minimum of a holonomic function $f(x)$ and the point $x$ which gives the minimum $f(x)$. In HGD, we use the iterative method together with a Pfaffian system. In this paper, we use the the Newton-Raphson iterative minimization method in which the renewal rule of the search point is
given by
$$
x_{k+1}=x_{k}-H^{-1}(x_{k})\nabla f(x_{k}),
$$
where $\displaystyle{ \nabla f(x_{k})=(\frac{\partial f(x_{k})}{\partial x_{1}},\ldots, \frac{\partial f(x_{k})}{\partial x_{k}})^{T}}$ and $H(x_{k})$ is the Hessian of $f(x)$ at $x=x_{k}$.

\subsection{Mathematical background}

HGD is based on the theory of the Gr\"{o}bner basis.
In the following, we refer to the relation of a numerical method and the Gr\"{o}bner basis.

Let $R$ be the differential ring written as
\begin{eqnarray}
R = {\bm C}[x_1,..., x_n] \langle \partial_1,..,\partial_n \rangle \nonumber
\end{eqnarray}
where ${\bm C}[x_1,...,x_n]$ are the rational coefficients of differential operators.
Suppose that  $I=\{\ell_i|i=1,...,p\} $ is a left ideal of $R$, $k[\bm{x}]$ is a field and $D \in k[\bm{x}]\langle \partial_1,..,\partial_n \rangle \in I$. If an arbitrary function $f$ satisfies $Df=0$ for all $D$, then $f$ is a solution of $I$. That is
\begin{eqnarray}\label{eq_h}
\ell_{i} f=0 \ \forall i
\end{eqnarray}
When $f$ satisfies Equation (\ref{eq_h}), $f$ is called {\it holonomic function}.

Let $s=[s_{1},\ldots,s_{t}]$,
with $s_{i}=q_{i}(x)\partial^{\alpha_{i}}$ be a standard basis in the quotient vector space $R/I$ which is a finite dimensional vector spaces.
Let $G$ be the Gr\"{o}bner basis of $I$. The rank of arbitrary differential operations can be reduced by normalization by $G$.
Assume that $\partial_{i} s_{j} \rightarrow_{G} \sum_{k}c^{i}_{jk}s_{k}$ holds. For a solution $f$ of $I$ put $F=(f,s_{2}f,\ldots,s_{t}f)^{T}$. Then, it holds that
\begin{prop}(see, e.g.,Nakayama et al.(2011))
\begin{eqnarray}
\frac{\partial F}{\partial x_{i}}=P_{i}F,\ \ i=1,...,n
\end{eqnarray}
where $P_{i}$ is a $t \times t$ matrix with $c^{i}_{jk}$ as  a $(j,k)$ element
\end{prop}
\begin{proof}
\begin{eqnarray}
\frac{\partial s_{j}f(x)}{\partial x_{i}} &=& (\partial_{i} \bullet s_{j})f(x)  \nonumber \\
                                          &=&  (\sum_{k}c^{i}_{jk}s_{jk})f(x) \ \  (mod \ \ I)  \nonumber \\
                                          &=&  [P_{i}F]_{j}, \ \ i=1,...,n,\ \ j=1...,t
\end{eqnarray}
This proves the assertion.
\end{proof}
The above differential equations are called {\it Pfaffian differential equations} or {\it Pfaffian system} of $I$. So we can calculate the gradient of $F$ by using Pfaffian differential equations. Then, $\nabla f(x_{k})$ and $H^{-1}(x_{k})$ are also given by Pfaffian differential equations. (see Hibi et al.(2012))
\begin{lemma}
Let $\sum_{j}^{t}a_{ij}s_{j}$ be the normal form of  $\partial_i=\partial / \partial x_i$ by $G$
and $\sum_{k}^{t}u_{ijk}s_{k}$ be the normal form of $\partial_i \partial_j$ by $G$. Then we have,
\begin{eqnarray}
\partial_{i} f(x_k) \rightarrow_{G} (\sum_{j}^{t}a_{ij}s_{j})f(x_k) =\sum_{j}^{t}a_{ij}F_j(x_k) = ((P_1F(x_k) )_1,...,(P_nF(x_k) )_1) \\
\partial_{i}\partial_{j}f(x_k) \rightarrow_{G} (\sum_{m}^{t}u_{ijk}s_{m})f(x_k)= \sum_{m}^{t}u_{ijm}F_m(x_k)=((\frac{\partial P_i}{\partial x_j}+P_iP_j)F(x_k))_1
\end{eqnarray}
where $(v)_1$ denotes the first entry of a vector $v$.
\end{lemma}

\subsection{Algorithm}

For HGD, we first give an ideal $I=\{\ell_i|i=1,...,p,\ \ell_{i} f=0 \ \forall i\} $ for holonomic function $f(x)$ and calculate the Gr\"{o}bner basis $G$ of $I$ and then the standard basis $S$ are given by $G$. The coefficient matrix $P_{i}$ for Pfaffian system is led by this standard basis, and  $H^{-1}(x_{k})$ and $\nabla f(x_{k})$ are calculated from $S$ by starting from a initial point $x_0$ through the Pfaffian equations. After these, we can compute automatically the optimum solution by a mixed use of  then Newton-Raphson method. The algorithm is given by below.
\newpage
\begin{algo}
\end{algo}\begin{itemize}\setlength{\itemindent}{1.5mm}
\item[step 1]Set $k=0$ and take an initial point $x_{0}$ and evaluate $F(x_0)=(f(x_0),s_{1}f(x_0),\ldots,s_{t}f(x_0))^{T}$.
\item[step 2]Evaluate $\nabla f(x_{k})$ and $-H^{-1}(x_{k})$ from $F$ and calculate the Newton direction, $d_k=-H^{-1}(x_{k})\nabla f(x_{k})$
\item[step 3]Update a search point by $x_{k+1}=x_{k}+d_{k}$.
\item[step 4]Evaluate $F(x_{k+1})$ by solving Pfaffian equations numerically.
\item[step 5]Set $k=k+1$ and calculate $F(x_{k+1})$ and  goes to step.2 and repeat until convergence.
\end{itemize}

\begin{remark} The key step of the above algorithm is step 4. We can not  evaluate $F(x_{k+1})$ by inputting $x_{k+1}$ in the function $f(x)$ since the HGD treats the case that $f(x)$ is difficult to calculate numerically.
Instead,  we only need calculate $f(x_0)$ and $F(x_0)$ numerically for a given initial value $x_0$.
\end{remark}

\section{Constrained holonomic gradient descent method}

Now, we propose the method in which we add constraint conditions to HGD and call it the constrained holonomic gradient descent method(CHGD).

\subsection{How to add the constraints}

For treating constraints  we use the penalty function and add it to objective function and make a new objective function and can treat it as the unconstrained optimization problem. We use HGD for evaluation of gradients and
Hessian  and use the exact penalty function method for constraints. The value of updating a search point can be obtained as the product of directional vector and step size. The step size $\alpha$ is chosen so that the following Armijo condition is satisfied. In fact we chose $\alpha$ such that
	\begin{eqnarray} \label{eq_s}
		P(x_k + \alpha x_k ; \rho) \leq P(x_k ; \rho) + \xi \alpha \{ P_l (x_k , \nabla x_k ; \rho) - P(x_k ; \rho) \},
	\end{eqnarray}
where $0 < \xi < 1$ and $P_l (x_k , \nabla x_k ; \rho)$ is the approximation of $P(x_k, \rho)$ given by.
	\begin{eqnarray}
		P_l (x_k , \nabla x_k ; \rho) &= &f(x_k) + \nabla f(x_k)^T \nabla x_k \nonumber \\
		& + & \rho \{  \sum_{i=1} ^m | g_i (x_k) + \nabla g_i (x_k) ^T \nabla x_k | \nonumber  \\
		& + &  \sum_{j=i} ^n max(0, h_j(x_k) + \nabla h_j (x_k) ^T \nabla x_k )\}
	\end{eqnarray}
The initial value of $\alpha$ is set $1$ and then $\alpha$ is made smaller iteratively until $\alpha$ satisfies Equation (\ref{eq_s}), or $\alpha=0$.

 In our algorithm, holonomic gradient descent plays a role to calculate the gradient vectors and then the penalty function plays a role to control the step size iteratively.

\section{Computational results}

We apply CHGD for MLE for von Mises distribution(vM). The process of applying for HGD is shown in Nakayama et al.(2011).
The density function of vM is given by  $f(\kappa, \mu)=e^{\kappa \cos (\mu - x)}/\int_0^{2\pi} e^{\kappa \cos( \mu - t)}dt$. The parameters of vM, $\kappa$ and $\mu$, show concentration and mean of angle data ${\bm x}$ respectively. We set the parameters for MLE $\theta_1=\kappa \cos \mu$ and $\theta_2 = \kappa \sin \mu $. Now we solve the constrained optimization problem given by.
\begin{eqnarray}\label{optvM}
(P)&min&\ \ L(\theta_1,\theta_2) = e^{-\overline{c}\theta_1-\overline{s}\theta_2}\int_{0}^{2\pi} e^{\theta_1 \cos t+\theta_2 \sin t} dt \nonumber \\
&s.t.&\  \theta_1 \leq \theta_2
\end{eqnarray}
Let  ${\bm x}$ be sample data. Let $n$ be sample size. Then, $\overline{c}=\frac{1}{n} \sum_i^n \cos x_i$ and $\overline{s}=\frac{1}{n} \sum_i^n \sin x_i$.

\subsection{Simulation}

In our simulation, we set the vM's parameter $(\kappa, \mu)=(5,\pi/4)$ of which the true value $(\theta_1, \theta_2)=(3.54, 3.54)$ and the initial value $(\theta_1, \theta_2)=(-2.0, 0.1)$. We tried the 2 patterns of constraints. Both of the case worked under the same condition except constraints.
In Figure 1, the constraint is $\theta_1 \leq \theta_2$. In Figure 2, the constraint is  $\theta_1^2+\theta_2^2 \leq 9$.
Figures 1,2 are the drawing of the trace of the search point.
\newline
\begin{figure}[htbp]
\begin{tabular}{cc}
\begin{minipage}{0.5\hsize}
\begin{center}
  \includegraphics[bb=70 0 533 539, scale=.4]{orbit4_2.png}
 \caption{the case of $ \theta_1 \leq \theta_2$}\end{center}
\end{minipage}
\begin{minipage}{0.5\hsize}
\begin{center}
  \includegraphics[bb=70 0 533 539, scale=.4]{orbit2_2.png}
 \caption{the case of $\theta_1^2+\theta_2^2 \leq 9$}
\end{center}
\end{minipage}
\end{tabular}
\end{figure}
\newline\newline
The result of simulation, the convergence point of HGD is $(\theta_1, \theta_2)=(3.63, 3.67)$.
In Figure 1, the convergence point of CHGD is $(\theta_1, \theta_2)=(0.13, -0.004)$. In Figure 2, the convergence point of CHGD is $(\theta_1, \theta_2)=(2.08,  1.94)$. In the CHGD, the search direction is almost same as the HGD, because the direction is decided by the HGD's algorithm.
While, the constraints play the role to judge the search point is within the feasible region or not and decide the step size.
\newpage

\subsection{Runtimes}

CHGD is the effective method for optimization with constraints.
However, whenever CHGD increases the cost of runtimes than HGD regardless of whether the solution is in the feasible region or not.
The following table shows the runtimes when the optimization solution is within the feasible region.
\newline
\begin{table}[htb]
  \begin{center}
  \caption{Comparing the runtimes\label{tb1}}
  \begin{tabular}{l||c|c|c} \hline
    & CPU TIME (sec) & \multicolumn{2}{c}{PATAMETERS $(\theta_1,\theta_2)$}   \\ \hline \hline
    HGD & 0.03698 &\ \ \ \ \ \ \ \ \ 2.120627\ \ \ \ \ \ \ \ \ \ & 2.120333 \\
    CHGD & 0.09834 & 2.120803 & 2.120629 \\
    NEWTON & 0.12598 & 2.124429 & 2.124855 \\\hline
    \multicolumn{4}{c}{\small We programmed by R and executed on Windows 7 64bit with RStudio Version 0.97.336}
  \end{tabular}
  \end{center}
\end{table}

In Table \ref{tb1}, all numbers are the means of 500 times trials. The optimization problem is Equation (\ref{optvM}). Sample data is drawn from the vM with $(\theta_1,\theta_2)=(2.12, 2.12)$. The third column of Table \ref{tb1} is the result with only Newton-Raphson method which optimize $f(x)$ directly, not use Pfaffian system. Thus, we see that HGD and CHGD is faster than Newton-Raphson method.

We see that the runtimes of CHGD is longer than HGD in general, where the both of solutions are almost the same value when the solution is inside  the feasible region. Sometimes the process finishes early by constraints, when the solution is outside the feasible region. Although, we need consider the cost of calculation of CHGD.