text:
\section{Mesh deformation techniques}\label{chap:methods}
The scope of this paper is to compare various techniques to construct ALE mappings. Put simply, an ALE mapping is nothing else but a deformation of the computational mesh in the fluid domain. In this section, we introduce several commonly used mesh deformation techniques (MDTs) as well as propose certain variations to them. All considered techniques achieve the same goal: provided a displacement of the FSI interface at a given time, they extend it into the fluid domain. Note that although we present these techniques in a 2D setting, one can readily apply them in 3D. When comparing different techniques, we largely pay our attention to the maximum amount of deformation a particular technique can handle. That is, how much the mesh can be deformed before the bijectivity condition (\ref{eq:bijectivity}) is violated. A secondary measure is of course the overall computational cost associated with computing ALE mappings using a given MDT.

After MDTs, we describe the Jacobian-based local stiffening \cite{stein2003mesh} which we use to augment each of the considered techniques. Finally, we discuss practical ways to check whether the bijectivity condition (\ref{eq:bijectivity}) is satisfied.

\subsection{Harmonic extension (HE/IHE)}
Probably the simplest way to extend displacement of the FSI interface into the fluid domain is by means of harmonic extension (HE) \cite{richter2010numerical,dorfel2011fluid,wu2014fully}. Given the interface displacement $\f{u}_s^i$ at time $t_i$, the ALE displacement $\f{u}_a^i$ is computed by solving Laplace's equation in the initial configuration of the fluid domain $\Omega_f^0$ for every displacement component:
\begin{align}
\Delta\f{u}_a^i &= \f{0}\text{ in }\Omega_f^0,\\
\f{u}_a^i &= \f{u}_s^i\text{ on }\Gamma^0,\\
\f{u}_a^i &= \f{0}\text{ on }\partial\Omega_f^0\setminus\Gamma^0.
\end{align}
The interface displacement $\f{u}_s^i$ serves as a Dirichlet boundary condition on the FSI interface $\Gamma^0$. At the rest of the boundary, the prescribed displacement is zero. Note that the method does not take into account information about the interface or ALE displacement from the previous time step.

The HE technique is the least computationally expensive method to construct ALE mappings. Let $N$ be the number of inner control points in the fluid domain. Since all displacement components satisfy the same equation, they can be computed by solving a single linear system with an $N\times N$ matrix and an $N\times d$ right-hand side. Here, $d$ is the problem dimension; in our case, $d=2$. Moreover, once the matrix is assembled, it can be reused for all time steps, which drastically reduces the computational cost associated with matrix assembly in IGA. At each time step, one only has to update the right-hand side to take the current FSI interface displacement into account. This can be performed efficiently by storing a Dirichlet elimination matrix when assembling the main matrix.

Despite its computational efficiency, the HE technique has serious disadvantages. First of all, it treats displacement components as completely independent variables and does not promote bijectivity of the ALE mapping in any way. Second, solutions of Laplace's equation in the vicinity of corners behave like $r^{\pi/\omega}$, where $r$ is the distance to the corner and $\omega$ is the corner angle. For reentrant corners, that is for $\omega>\pi$, solutions do not belong to $H^1(\Omega_f^0)$ since their derivatives tend to infinity. As a consequence, the corresponding ALE mapping may lose its bijectivity. Due to these two problems, the HE technique usually is only able to handle rather small deformations.

We propose a slight improvement to the HE technique achieved by turning it into an incremental algorithm. Assume that the ALE displacement $\f{u}_a^i$ at time $t_i$ is known. We can use it to deform the initial configuration of the fluid domain $\Omega_f^0$ and obtain the deformed configuration $\Omega_f^i = \Omega_f^0 + \f{u}_a^i$ as equation (\ref{eq:deformDomain}) describes. We can then compute an ALE displacement increment $\delta\f{u}_a^{i+1}$ by solving Laplace's equation in the deformed configuration:
\begin{align}
\label{eq:ihe1}\Delta\delta\f{u}_a^{i+1} &= \f{0}\text{ in }\Omega_f^i,\\
\delta\f{u}_a^{i+1} &= \f{u}_s^{i+1}-\f{u}_s^i\text{ on }\Gamma^i,\\
\label{eq:ihe2}\delta\f{u}_a^{i+1} &= \f{0} \text{ on }\partial\Omega_f^i\setminus\Gamma^i.
\end{align}
And finally, we define the ALE displacement $\f{u}_a^{i+1}$ at time $t_{i+1}$ as $\f{u}_a^i+\delta\f{u}_a^{i+1}$. Note that the resulting incremental harmonic extension (IHE) technique is not equivalent to the HE technique since ALE increments are computed in deformed configurations of the fluid domain.

One advantage of the IHE technique is that it uses the ALE displacement from the previous time step. Therefore, IHE can be expected to perform slightly better than the HE technique, meaning that it can handle larger deformations. On the other hand, each step of the IHE technique is formulated in a different configuration of the fluid domain than the previous one. As a result, the technique requires matrix assembly at every time step, which makes it more computationally expensive than HE.

In what follows, we apply the same ideas to other MDTs and consider both their non-incremental and incremental versions, which often share the same advantages and disadvantages as the HE and IHE techniques.

\subsection{Bi-harmonic extension (BE/IBE)}
The HE technique is weak when it comes to large deformations. To overcome this problem, one can search for the ALE displacement as a solution to the bi-harmonic equation:
\begin{align}
\label{eq:biharmeq}\Delta^2\f{u}_a &= \f{0}\text{ in }\Omega_f.
\end{align}
Solutions of the bi-harmonic equation are known to be more regular in comparison to Laplace's equation and do not have problems at reentrant corners \cite{richter2010numerical}.

However, this bi-harmonic extension (BE) technique is often dismissed as too computationally expensive. Indeed, in order to solve the bi-harmonic equation, one has two options: either use $C^1$-conforming elements, which in IGA requires $G^1$-continuity between patches \cite{collin2016analysis,birner2018approximation}; or use mixed elements with an auxiliary variable $\f{q}$ to replace the bi-harmonic equation with two Laplace's equations \cite{boffi2013mixed}:
\begin{equation}
\Delta\f{u}_a = \Delta\f{q},\sg\Delta\f{q}= \f{0}\text{ in }\Omega_f.
\end{equation}
In this work, we consider only the latter option since it is easier to implement for multi-patch geometries.

In our interpretation, the BE technique has the following form: given the interface displacement $\f{u}_s^i$ at time $t_i$, the ALE displacement $\f{u}_a^i$ is computed by solving the following linear system in the initial configuration of the fluid domain $\Omega_f^0$:
\begin{align}
\label{eq:BE1}\Delta\f{u}_a^i &= \Delta\f{q},\sg \Delta \f{q} = \f{0}\text{ in }\Omega_f^0,\\
\f{u}_a^i &= \f{u}_s^i \text{ on }\Gamma^0,\\
\f{u}_a^i &= \f{0}\text{ on }\partial\Omega_f^0\setminus\Gamma^0,\\
\label{eq:BE2}\nabla\f{q}\cdot\f{n} &= \f{0} \text{ on }\partial\Omega_f^0.
\end{align}

The BE technique shares many similarities with HE which make both techniques very efficient: it does not use information from previous time steps; the linear system has to be assembled only once; the multiple-right-hand-sides approach can be used to compute all displacement components at once. However, the resulting linear system is twice the size of the HE linear system: the matrix is of size $2N\times2N$, and the right-hand size is of size $2N\times d$. Moreover, the linear system has a saddle-point structure, so specialized linear solvers are necessary to solve it efficiently.

Just like with the HE technique, we propose an incremental variation of the bi-harmonic extension (IBE). An increment $\delta\f{u}_a^{i+1}$ is computed by solving the following system in the deformed configuration $\Omega_f^i$:
\begin{align}
\Delta(\delta\f{u}_a^{i+1}) &= \Delta\f{q},\sg \Delta \f{q} = \f{0}\text{ in }\Omega_f^i,\\
\delta\f{u}_a^{i+1} &= \f{u}_s^{i+1}-\f{u}_s^i\text{ on }\Gamma^i,\\
\delta\f{u}_a^{i+1} &= \f{0}\text{ on }\partial\Omega_f^i\setminus\Gamma^i,\\
\nabla\f{q}\cdot\f{n} &= \f{0} \text{ on }\partial\Omega_f^i.
\end{align}
After that, the ALE displacement $\f{u}_a^{i+1}$ at time $t_{i+1}$ is defined as $\f{u}_a^i+\delta\f{u}_a^{i+1}$.

The IBE technique requires matrix assembly at each time step but can potentially handle larger mesh deformations than the BE technique.

\subsection{Linear elasticity (LE/ILE)}
The next MDT we consider is based on the linear elasticity theory. It is widely used in FSI applications and belongs to the state-of-the-art in the field \cite{bazilevs2013,stein2003mesh,tezduyar2007modelling}. The core idea is to treat the fluid domain as an elastic body and to construct ALE displacement as a solution to the equations of linear elasticity:
\begin{equation}
\dv\pmb{\sigma}_a(\f{u}_a) = \f{0}\text{ in }\Omega_f.
\end{equation}
Here, $\pmb{\sigma}_a$ is the Cauchy stress tensor related to the linearized strain tensor $\pmb{\varepsilon}_a = (\nabla\f{u}_a+\nabla\f{u}_a^T)/2$ by the Hooke's law:
\begin{equation}
\pmb{\sigma}_a = \lambda_a\tr(\pmb{\varepsilon}_a)\f{I} + 2\mu_a\pmb{\varepsilon}_a.
\end{equation}
The Lam\'e parameters $\lambda_a$ and $\mu_a$ depend on Young's modulus $E_a$ and Poisson's ratio $\nu_a$. Since we do not apply volumetric or surface force to the fluid domain, Young's modulus does not affect the resulting ALE displacement. On the other hand, Poisson's ratio is important because it regulates resistance of the fluid mesh to volumetric changes. A too high value (close to 0.5) would result in an almost incompressible behavior, which could lead to excessive distortion of the mesh elements and numerical instabilities. In contrast to that, a too low value (close to 0 or even negative) can reduce resistance of the fluid mesh to bijectivity violation. Therefore, we recommend choosing a value between 0.3 and 0.45.

Unlike the techniques based on harmonic and bi-harmonic extension, the linear elasticity MDT is best known in its incremental version. That is, given the ALE displacement $\f{u}_a^i$ at time $t_i$, an increment $\delta\f{u}_a^{i+1}$ is computed by solving the linear elasticity equations in the deformed configuration of the fluid domain $\Omega_f^i$:
\begin{align}
\label{eq:ILE1}\dv\pmb{\sigma}_a(\delta\f{u}_a^{i+1}) &= \f{0}\text{ in }\Omega_f^i,\\
\delta\f{u}_a^{i+1} &= \f{u}_s^{i+1}-\f{u}_s^i\text{ on }\Gamma^i,\\
\label{eq:ILE2}\delta\f{u}_a^{i+1} &= \f{0}\text{ on }\partial\Omega_f^i\setminus\Gamma^i.
\end{align}
After that, the ALE displacement $\f{u}_a^{i+1}$ at time $t_{i+1}$ is defined as $\f{u}_a^i+\delta\f{u}_a^{i+1}$. We refer to this technique as incremental linear elasticity (ILE). The ILE technique is known for its robustness and an ability to withstand large mesh deformations. However, little to no research has been conducted to explain its superior behavior.

With respect to computational cost, the ILE technique involves solving a linear system with a matrix of size $dN\times dN$ and a right-hand size of size $dN\times1$. The linear system has to be reassembled at each time step. Note that the size of the linear system depends on a dimension of the problem. Therefore, it scales worse from 2D to 3D than the IHE and IBE techniques.

For the sake of completeness, let us also study a non-incremental version of the ILE technique. The ALE displacement $\f{u}_a^i$ at time $t_i$ is computed by solving the equations of linear elasticity in the initial configuration of the fluid domain $\Omega_f^0$:
\begin{align}
\dv\pmb{\sigma}_a(\f{u}_a^i) &= \f{0}\text{ in }\Omega_f^0,\\
\f{u}_a^i &= \f{u}_s^i\text{ on }\Gamma^0,\\
\f{u}_a^i &= \f{0} \text{ on }\partial\Omega_f^0\setminus\Gamma^0.
\end{align}
We call this the linear elasticity (LE) technique. Although one can only expect it to perform well for small deformations, it is rather computationally inexpensive. Similarly to the HE and BE techniques, the LE technique requires matrix assembly only once and lets one reuse the matrix for every time step.

\subsection{Nonlinear elasticity (TINE)}
The last MDT we present in this paper is based on equations of nonlinear elasticity. The idea is to construct the ALE displacement at each time step as an approximate solution to the local balance equations of linear momentum
\begin{equation}\label{eq:nonlinel_tine}
\dv \f{P}_a(\f{u}_a^i) = \f{0}\text{ in }\Omega_f^0,
\end{equation}
where $\f{P}_a = \f{F}_a\f{S}_a$. To ensure bijectivity of the ALE mapping, we use a logarithmic variation of the neo-Hookean material law
\begin{equation}\label{eq:neohooke}
\f{S}_a = \lambda_a\ln J_a\f{C}_a^{-1} + \mu_a(\f{I}-\f{C}_a^{-1}),
\end{equation}
where $\f{C}_a = \f{F}_a^T\f{F}_a$ is the right Cauchy-Green strain tensor. Similarly to the LE and ILE techniques, the Lam\'e parameters $\lambda_a$ and $\mu_a$ can be computed from Young's modulus $E_a$ and Poisson's ratio $\nu_a$, of which only Poisson's ratio affects the solution of equations (\ref{eq:nonlinel_tine}).

Due to the term $\ln J_a$  in the neo-Hookean law (\ref{eq:neohooke}), any solution of equations (\ref{eq:nonlinel_tine}) satisfies the bijectivity condition (\ref{eq:bijectivity}). This fact makes the MDT based on equations (\ref{eq:nonlinel_tine}) unique since it explicitly enforces the bijectivity condition. Unfortunately, equations (\ref{eq:nonlinel_tine}) are nonlinear, and an attempt to fully solve them at each time step would make the MDT prohibitively expensive. However, since the ALE mapping should possess certain regularity in time, it is possible to use a solution of equations (\ref{eq:nonlinel_tine}) at time $t_i$ to efficiently construct an approximate solution at time $t_{i+1}$ \cite{shamanskiy2019mesh}. We refer to this technique as tangential incremental nonlinear elasticity (TINE). The TINE technique can be seen as pseudo time-stepping or an example of the continuation methods for nonlinear problems \cite{deuflhard2011newton}.

Let us look under the hood of TINE. It is based on a Newton-like linearization of equations (\ref{eq:nonlinel_tine}). To define it, we need to transform equations (\ref{eq:nonlinel_tine}) into a weak form, also known as variation formulation. To that end, let us define a solution space $\m{V} = (H^1(\Omega_f^0))^d$ and a test space $\m{V}_0 = \{\f{w}\in (H^1(\Omega_f^0))^d\;|\;\f{w}=\f{0}\text{ on }\partial\Omega_f^0\}$. We can then write the weak form of equations (\ref{eq:nonlinel_tine}) as
\begin{align}
\nonumber&\text{find } \f{u}_a\in\m{V}\text{ such that }\forall\f{w}\in\m{V}_0\\
&R(\f{u}_a,\f{w}) = \int\displaylimits_\Omega\f{S}_a:\delta\f{E}_a[\f{w}]d\f{x} = 0.
\end{align}
Here, $\delta\f{E}_a[\f{w}] = \frac{1}{2}\Big(\f{F}_a^T\nabla\f{w} + \nabla\f{w}^T\f{F}_a\Big)$ is the variation of the Green-Lagrange strain tensor. The Taylor expansion at point $(\f{u}_a,\f{w})$ with an increment $\delta\f{u}_a$ yields
\begin{align}
\nonumber R(\f{u}_a+\delta\f{u}_a,\f{w}) = &R(\f{u}_a,\f{w}) +\\ &DR(\f{u}_a,\f{w})\cdot\delta\f{u}_a + o(||\delta\f{u}_a||),
\end{align}
where $DR(\f{u}_a,\f{w})\cdot\delta\f{u}_a$ is a directional derivative. We refer to \cite{wriggers2008nonlinear,shamanskiy2019mesh} for details on computing $DR(\f{u}_a,\f{w})\cdot\delta\f{u}_a$.

The idea of the TINE technique is to use one Newton-like step
\begin{align}
\nonumber&\text{find } \delta\f{u}_a\in\m{V}\text{ such that }\forall\f{w}\in\m{V}_0\\
&DR(\f{u}_a,\f{w})\cdot\delta\f{u}_a = -R(\f{u}_a,\f{w})
\end{align}
per time step to compute an ALE increment and update the ALE displacement. Concretely, given the ALE displacement $\f{u}_a^i$ at time $t_i$, we find an ALE increment $\delta\f{u}_a^{i+1}$ as a solution of the linear problem
\begin{align}
DR(\f{u}_a^i,\f{w})\cdot\delta\f{u}_a^{i+1} &= -R(\f{u}_a^i,\f{w}) \sg\forall\f{w}\in\m{V}_0,\\
\delta\f{u}_a^{i+1} &= \f{u}_s^{i+1}-\f{u}_s^i\text{ on }\Gamma^0,\\
\delta\f{u}_a^{i+1} &= \f{0} \text{ on }\partial\Omega_f^0\setminus\Gamma^0.
\end{align}
After that, we define the ALE displacement $\f{u}_a^{i+1}$ at time $t_{i+1}$ as $\f{u}_a^i+\delta\f{u}_a^{i+1}$.

It is natural to compare the TINE and ILE techniques which are very similar at first glance. Both are incremental techniques based on the elasticity theory; both require solution of a linear system with a matrix of size $dN\times dN$ and a right-hand side of size $dN\times1$; both require matrix assembly at each time step. In general, one can expect both techniques to be roughly equal in computational cost. Unlike ILE, however, the TINE technique explicitly enforces the bijectivity condition (\ref{eq:bijectivity}). Moreover, the TINE technique is based in the initial configuration of the fluid domain. As we show in Sections \ref{chap:testALE} and \ref{chap:testFSI}, this last observation results in crucial differences in behavior of the ILE and TINE techniques when it comes to the accumulated distortion effect.

\subsection{Local stiffening}
Most of the fluid mesh deformation happens along the FSI interface, where the structure displacement is applied as a Dirichlet boundary conditions to the ALE displacement. On the other hand, mesh elements in the vicinity of the stationary boundary of the fluid domain undergo almost no deformation. Therefore, their contribution into processing of the applied interface displacement is negligible. If the deformation could be redistributed away from the FSI interface towards the stationary boundary, the mesh could undergo larger deformations without becoming invalid. This is the idea behind local stiffening, which locally changes the way different elements react to the deformation.

Let $\f{G}:[0,1]^d\to\Omega$ be a parametrization of the computational domain $\Omega$. Imagine that we have to compute integrals corresponding to matrix entries of the discretized linear system. One of the simplest ways to implement local stiffening is to drop the Jacobian determinant $\det\nabla\f{G}$ when transforming the integrals from domain $\Omega$ to parametric domain $[0,1]^d$:
\begin{equation}\label{eq:simplelocstiff}
\int\displaylimits_\Omega(\cdots)d\f{x} = \int\displaylimits_{[0,1]^d}(\cdots)\det\nabla\f{G}d\pmb{\xi}\to \int\displaylimits_{[0,1]^d}(\cdots)d\pmb{\xi}.
\end{equation}
This local stiffening method was first proposed in \cite{tezduyar1993parallel}. For elasticity problems, (\ref{eq:simplelocstiff}) can be interpreted as a local change of Young's modulus
\begin{equation}\label{eq:fsi:stiffE}
E\to\frac{E}{\det\nabla_{\pmb{\xi}}\f{G}},
\end{equation}
which makes elements with small values of $\det\nabla\f{G}$ stiffer and elements with large values softer. Therefore, the former elements undergo less deformation and are less likely to become invalid.

A more advanced local stiffening method introduced in \cite{stein2003mesh} does not simply drop the Jacobian determinant but changes the degree with which it enters the integrals:
\begin{align}\label{eq:localstiff}
\nonumber\int\displaylimits_\Omega(\cdots)d\f{x} = &\int\displaylimits_{[0,1]^d}(\cdots)\det\nabla\f{G}d\pmb{\xi}\to\\
&\int\displaylimits_{[0,1]^d}(\cdots)(\operatorname{det}\nabla\f{G})^{1-\chi}d\pmb{\xi}.
\end{align}
Here, $\chi\geqslant0$ is called the stiffening degree. The higher the stiffening degree is, the more local stiffening is achieved. $\chi=0$ corresponds to no local stiffening, and $\chi=1$ corresponds to Jacobian dropping (\ref{eq:simplelocstiff}). Too high stiffening degrees, however, may result in excessive mesh distortion. We refer to this method as the Jacobian-based local stiffening.

The Jacobian-based local stiffening acts differently depending on whether the mesh deformation method is formulated in the initial or in the deformed configuration of the fluid domain. Namely, if the integrals for matrix entries are computed in the initial configuration $\Omega_f^0$, the local stiffening is based only on the initial parametrization $\f{G}_f^0$. However, in the case of the deformed configuration $\Omega_f^i$, the local stiffening takes into account already applied deformation since the parametrization $\f{G}_f^i$ of the deformed configuration is defined as $(\f{I}+\f{u}_a^i)\circ\f{G}_f^0$, see equation (\ref{eq:deformDomain}). This effect has both advantages and disadvantages. From one point of view, if a particular mesh element becomes ill-shaped after the deformation, its value of $\det\nabla\f{G}_f^i$ decreases. As a result, this element receives more stiffening, which prevents it from becoming even more ill-shaped or invalid. On the other hand, in case of the ILE technique, this deformation-aware local stiffening essentially makes material properties of the mesh deformation-dependent, which can cause irreversible plastic deformation. For other MDTs based in the deformed configuration, namely IHE and IBE, the effect is similar. As we show in Sections \ref{chap:testALE} and \ref{chap:testFSI}, this irreversible deformation accumulates over time and can significantly affect results of FSI simulations. We refer to this effect as accumulated distortion.

Regardless if the local stiffening is deformation-aware or not, the initial parametrization $\f{G}_f^0$ of the fluid domain provides a major contribution to how much stiffening each element receives. Let us consider the deforming part of the fluid domain, see Figure \ref{fig:localstiff}. The top and bottom patches are perfect parallelograms, and $\det\nabla\f{G}_f^0$ is constant. Therefore, elements of these patches receive no local stiffening with respect to each other. However, the right patch has a distinct tapered left side, where $\det\nabla\f{G}_f^0$ becomes very small in comparison to the surrounding elements of the top and bottom patches. As a result, element on the left side of the right patch become much stiffer and maintain their shape. Consequently, angles of all three mesh patches that are adjacent to the beam right end do not change much during mesh deformation. In particular, they do not exceed $\pi$, which would lead to $\det\nabla\f{G}_f^i$ becoming negative, which means that the bijectivity condition (\ref{eq:bijectivity}) is not violated.

\begin{figure}[H]
	\centering
	\includegraphics[height=4cm]{images/fig3.png}
	\caption{Local stiffening potential of the deforming fluid mesh part.}
	\label{fig:localstiff}
\end{figure}

Note that local stiffening may be harder to achieve in IGA than in classical FEM. We have seen an example of a parallelogram with a natural uniform tensor-product NURBS parametrization. In that case, $\det\nabla\f{G}$ is constant throughout the geometry, and the Jacobian-based local stiffening would have no effect. With a finite element discretization of a parallelogram, it would be possible to place smaller elements where necessary. Since mappings to a reference element are independent for each FEM element, the Jacobian-based local stiffening would have an effect.

In IGA, it is possible to locally control the element size distribution with local refinement, for example by means of THB-splines \cite{vuong2011hierarchical}. However, local refinement does not change the underlying geometry parametrization and values of the Jacobian determinant. Therefore, in order to apply the Jacobian-based local stiffening in IGA, one has two options: either use non-uniform parametrizations with carefully designed knot vectors to artificially construct domain regions with smaller values of the Jacobian determinant; or design patch geometries in a way that naturally defines such domain regions. An example in Figure \ref{fig:localstiff} belongs to the latter approach.

\subsection{Bijectivity check}
Let us briefly discuss ways to check the bijectivity condition (\ref{eq:bijectivity}) in practice. A solution which takes the NURBS nature of the ALE displacement $\f{u}_a$ into account is to express $J_a$ as a NURBS function \cite{Gravesen2014}. If all coefficients in a NURBS representation of $J_a$ are positive, then the displacement $\f{u}_a$ satisfies the bijectivity condition. Unfortunately, this condition is only sufficient and not a necessary one. Therefore, it may often lead to false detection of bijectivity violation. In practice, we resort to a less elegant solution of sampling $J_a$ at the Gaussian quadrature points associated with the NURBS basis of $\f{u}_a$.

Note that whichever method is chosen, it introduces a certain computational overhead to construction of ALE mappings. Nevertheless, we recommend doing some bijectivity check at every time step, or at least with regular intervals. An ALE mapping that does not satisfy the bijectivity condition (\ref{eq:bijectivity}) quickly makes all ensuing computation results meaningless.


TITLE: Reciprocity map in global class field theory
QUESTION [2 upvotes]: Say $K$ is a number field, $v$ is a finite place of $K$, $K_v$ the $v$-adic completion of $K$.
We have the local Artin map for every finite $v$:
$$\rho_v : K_v^{\times}\to\text{Gal}(K_v^{\rm ab}/K_v)$$
sending $\varpi_v$, a generator of the copy of $\mathbf{Z}$ in $K_v^{\times}$, to the Frobenius element in $\text{Gal}(K_v^{\rm unr}/K_v)$, and sending $\mathcal{O}_v^{\times}$ isomorphically onto $\text{Gal}(K_v^{\rm LT}/K_v)$, where $K_v^{\rm LT}$ is the Lubin-Tate extension of $K_v$, with $K_v^{\rm LT}\cdot K_v^{\rm unr} = K_v^{\rm ab}$.
The completion map $K\to K_v$ induces an injective map $$\alpha_v : \text{Gal}(K_v^{\rm ab}/K_v)\to\text{Gal}(K^{\rm ab}/K)$$
Let's call $c_v$ the map $K^{\times}\to K_v^{\times}$ induced by completion.
Is there an explicit description for:
$$R_v := \alpha_v\circ\rho_v\circ c_v : K^{\times}\to\text{Gal}(K^{\rm ab}/K)\ ?$$

REPLY [1 votes]: In view of the so called "global-local principle", it does not seem very natural to consider only one single place $v$ of $K$. Usually, one fixes a rational prime $p$ and the set $S$ of all places of $K$ above $p$, and one introduces the maximal abelian extension $K_S$ of $K$ which is unramified outside $S$, as well as the subfield $H_S$ of $K_S$ in which all the places of $S$ are completely decomposed.  Note that $H_S$ is a modified Hilbert class field, with $Gal(H_S /K) \cong Cl_S (K)$, the quotient of the class group of $K$ modulo the subgroup generated by the classes of the primes in $S$. Then CFT (in idelic terms) gives a natural exact sequence $1 \to \bar E_S \to \prod_S Gal({K_v}^{ab}/K_v) \to Gal(K_S /K) \to Gal(H_S /K) \to 1$. The leftmost term is the closure of the group of $S$-units of $K$ embedded diagonally in the product of the profinite completions $\bar{K_v}^*$ of the groups ${K_v}^*, v\in S$, each $\bar{K_v}^*$ being $\cong Gal({K_v}^{ab}/K_v)$ via $\rho_v$. If we wave hands and say that we "know" the two extreme terms, then we "know"  what you call the maps $R_v$ provided we know the local $\rho_v$'s ("local reciprocity laws" of Hasse, Brückner, Kato, etc.). Note that your maps $\alpha_v$ cannot be injective in this setting.


text:
\begin{document}
\theoremstyle{plain}
\newtheorem{thm}{Theorem}[section]
\newtheorem{lem}[thm]{Lemma}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{prop}[thm]{Proposition}

\title{Phase transitions in non-linear urns with interacting types}
\author{Marcelo Costa \\ Universidad de Buenos Aires \thanks{\href{mailto:mrocha@dm.uba.ar}{mrocha@dm.uba.ar}} \and Jonathan Jordan \\ University of Sheffield \thanks{\href{mailto:jonathan.jordan@sheffield.ac.uk}{jonathan.jordan@sheffield.ac.uk}}}
\maketitle

\begin{abstract}We investigate reinforced non-linear urns with interacting types, and show that where there are three interacting types there are phenomena which do not occur with two types. In a model with three types where the interactions between the types are symmetric, we show the existence of a double phase transition with three phases: as well as a phase with an almost sure limit where each of the three colours is equally represented and a phase with almost sure convergence to an asymmetric limit, which both occur with two types, there is also an intermediate phase where both symmetric and asymmetric limits are possible.  In a model with anti-symmetric interactions between the types, we show the existence of a phase where the proportions of the three colours cycle and do not converge to a limit, alongside a phase where the proportions of the three colours can converge to a limit where each of the three is equally represented.\end{abstract}

\section{Introduction and definitions}

In general terms, an urn model is a system containing a number of particles of different types
(often regarded as balls of different colours, for ease of visualisation).
At each time step, a set of particles is sampled from the system, whose contents are then altered depending on the sample which was drawn.
Pemantle \cite{pemantlesurvey} surveys several ways to approach this model framework.

This paper is limited to models with a single urn from which a single ball is drawn, its colour is noted and it is then returned to the urn along with one new ball of that same colour.
In addition, we introduce a graph-based interaction according to which the probability of choosing a ball of a given colour is reinforced not only by its own proportion in the urn, but also by the proportions of balls of other colours.
Therefore, the interaction arises among balls of different colours, as opposed to the so-called interacting urn models consisting of  systems of multiple urns (e.g. Bena\"{i}m et al \cite{BBCL15}, Launay and Limic \cite {LL12}) in which different urns (each containing balls of different colours) interact with one another.  Our model is also different from the graph-based competition described in van der Hofstad et al \cite{van2016strongly}, where the colours correspond to edges of the graph, which compete with, as opposed to being reinforced by, other edges incident on the same vertices.

We now formally define our model.  Consider an urn containing balls of $d$ colours.
The vector $x(n) = (x_1(n),\ldots, x_d(n)) \in \N^d$ denotes the number of balls of each colour at time $n=0,1,2,...$.
The strength of the reinforcement is given by a positive real number $\beta>0$ and we denote by $x^{\beta}(n)$ the coordinate-wise $\beta$ power of the column vector $x(n)$.
The interaction is defined as follows.
Given a non-negative matrix $A = (a_{ij})_{i,j=1}^d$, define the column vector

\begin{equation}
  \label{eq:potential}
  u(n) := A x^{\beta}(n)
\end{equation}
Let $\salj_n$ be the $\sigma$-algebra generated by the $x(m)$ for $0\leq m\leq n$ and write $u_i(n)$ for the $i$-th component of $u(n)$.
The transition probabilities are then
\begin{equation}
  \label{eq:transition}
\pr(x(n+1)-x(n)=\mathbf{e}_i \> | \> \salj_n)=\frac{u_i(n)}{\sum_{j=1}^{N} u_j(n)}, \quad i=1,\ldots,d,
\end{equation}
where $\mathbf{e}_i$ is the unit vector in direction $i$.  That is, one ball is added to the urn at each time step, and the right hand side of \eqref{eq:transition} gives the probability that it is of colour $i$.

Now, let $n_0$ be the initial number of balls so that at time $n$ the urn contains $n+n_0$ balls.
Then, the proportion of balls of each colour is a process in the $(d-1)$-dimensional simplex $\Delta^{d-1}$ given by the vector
\begin{equation}
  \label{eq:proportion}
\bar{x}(n)=x(n)/(n+n_0).
 \end{equation}
Assuming that $A$ is a multiple of the identity matrix (therefore, $\bar{x}(n)$ having no interaction), it is well-known (see Oliveira \cite{RO07}) that the process $\bar{x}$ undergoes a phase transition as follows.
For $\beta < 1$, the process converges almost surely to the `centre' of the simplex, that is, the asymptotic proportion of balls of each colour are all the same.
For $\beta = 1$, commonly referred to as the P\'olya urn model, the process converges almost surely to a non-trivial random variable supported in the interior of the simplex.
For $\beta > 1$, the process converges almost surely to one of the corners of the simplex.
In this case that a single type dominates was proved by Khanin and Khanin \cite{khanin2001} following on from the two-type case which can be covered using Rubin's Theorem in Davis \cite{davis_rrw}.

For a two-colour urn model $d=2$ and symmetric interaction,
\[   A =
 \left(
 \begin{matrix}
   1 & a \\
   a & 1
 \end{matrix}
\right), \quad a>0, \]
it was proved by the first author in Theorem 2.2.1, \cite{MCthesis}, that there was a phase transition as follows.
\begin{align}
  & (i) \quad  \text { if } \quad  \left(\tfrac{1-a}{1+a}\right) \beta \leq 1, \quad \text{then} \quad  \bar{x}(n) \rightarrow (\tfrac{1}{2},\tfrac{1}{2}) \quad a.s. \nonumber\\
  & (ii) \quad \text {if } \quad  \left(\tfrac{1-a}{1+a}\right) \beta  > 1,   \quad \text{then} \quad  \bar{x}(n) \rightarrow  \Psi \quad a.s., \nonumber
\end{align}
where $\Psi$ is a random vector supported on $\left\{ \left( \frac{1}{1+r}, \frac{r}{1+r} \right), \left( \frac{r}{1+r}, \frac{1}{1+r}\right)\right\}$
and $r:= r(a, \beta)$ is the unique root in $(0,1)$ of $\mathcal{P}_{a,\beta}(z) = az^{\beta+1}- z^{\beta}+ z -a = 0.$
In case (\emph{ii}),  $\pr[\bar{x}(n) \rightarrow (\frac{1}{2},\frac{1}{2})] = 0.$
Note that for $\beta=1$, the process $(u(n))_{n \geq 0}$ is a Friedman's urn model and statement (\emph{i}) yields $u(n)/(u_1(n)+u_2(n)) \rightarrow (\frac{1}{2}, \frac{1}{2})$ $a.s.$ as expected.

In this paper, we follow up on the results for $d=2$ in \cite{MCthesis}, with the aim to generalise from $d=2$ to larger values of $d$ and to see whether more types of behaviour emerge when this is done.  We show that this is indeed the case when $d=3$, where we consider two particular choices of $A$.

First of all, we consider a choice of $A$ with a symmetric interaction of the same strength $a$ for each pair of colours.    The following theorem shows that in this system there are three phases, as opposed to two when $d=2$; there are phases where there is almost sure convergence to a symmetric limit and where there is almost sure convergence to one of a number of asymmetric limits, which are analogues of the phases when $d=2$, but there is also an intermediate phase where both symmetric and asymmetric limits are possible.

\renewcommand{\labelenumi}{(\roman{enumi})}
\begin{thm}\label{sym3main}
Let $A$ be the matrix
\begin{equation}
  \begin{pmatrix} 1 & a & a \\ a & 1 & a \\ a & a & 1 \end{pmatrix}, \quad a > 0.
\end{equation}\begin{enumerate}
\item Fix $a<1$.  Then there exists $\beta_1(a)$ satisfying $1<\beta_1(a)<\frac{1+2a}{1-a}$, with $\beta_1(a)$ an increasing function of $a$ satisfying $\beta_1(a)\to\infty$ as $a\to 1$, and we have the following three phases.
\begin{enumerate}
\item \label{symlim} \emph{Symmetric limit almost surely.} If $\beta<\beta_1(a)$, then almost surely $\bar{x}(n) \to \left(\frac13,\frac13,\frac13\right)$.
\item \label{bothlim} \emph{Symmetric or asymmetric limit.} If $\beta_1(a)<\beta<\frac{1+2a}{1-a}$ then there exists $r_2>1$ such that almost surely $\bar{x}(n)$ converges to one of the four points in $\Delta^2$ given by $\left(\frac13,\frac13,\frac13\right)$, $\left(\frac{r_{2}}{2+r_{2}},\frac{1}{2+r_{2}},\frac{1}{2+r_{2}}\right)$,  $\left(\frac{1}{2+r_{2}},\frac{r_{2}}{2+r_{2}},\frac{1}{2+r_{2}}\right)$ and  $\left(\frac{1}{2+r_{2}},\frac{1}{2+r_{2}},\frac{r_{2}}{2+r_{2}}\right)$.  All of these points have positive probability of being limits.
\item \label{asymlim} \emph{Asymmetric limit almost surely.} If $\beta>\frac{1+2a}{1-a}$ then there exists $r_{+}>1$ such that almost surely $\bar{x}(n)$ converges to one of the three points in $\Delta^2$ given by $\left(\frac{r_{+}}{2+r_{+}},\frac{1}{2+r_{+}},\frac{1}{2+r_{+}}\right)$,  $\left(\frac{1}{2+r_{+}},\frac{r_{+}}{2+r_{+}},\frac{1}{2+r_{+}}\right)$ and  $\left(\frac{1}{2+r_{+}},\frac{1}{2+r_{+}},\frac{r_{+}}{2+r_{+}}\right)$.
\end{enumerate}
\item Fix $a\geq 1$.  Then almost surely $\bar{x}(n) \to \left(\frac13,\frac13,\frac13\right)$.
\end{enumerate}
\end{thm}

Theorem \ref{sym3main} presents the results in terms of phase transitions in $\beta$ with $a$ fixed.  However, because both $\beta_1(a)$ and $\frac{1+2a}{1-a}$ are increasing functions of $a$ which converge to $1$ as $a\to 0$ and to $\infty$ as $a\to 1$, it is also possible to see them as phase transitions in $a$ with $\beta>1$ fixed: if $a<\frac{\beta-1}{2+\beta}$ then we will be in case (c), if $\frac{\beta-1}{2+\beta}<a<\beta_1^{-1}(\beta)$ then we will be in case (b), and if $a>\beta_1^{-1}(\beta)$ we will be in case (a).

We also consider a system where each colour is reinforced by itself and by one other, in a cyclic way.  For this system, the following theorem shows the existence of a phase transition between a phase with convergence with positive probability to a symmetric limit and a phase where there is no convergence to a limit and there is cycling behaviour.

\begin{thm}\label{cyclicmain}Let $A$ be the matrix
$$\begin{pmatrix} 1 & 1 & 0 \\ 0 & 1 & 1 \\ 1 & 0 & 1 \end{pmatrix}.$$
\begin{itemize} \item When $\beta<4$ there is positive probability that $\bar{X}(n)\to (1/3,1/3,1/3)$. \item When $\beta>4$, almost surely $\bar{X}(n)$ fails to converge, and the limit set is either a periodic orbit or a connected union of periodic orbits.\end{itemize}\end{thm}

In Section \ref{stochapprox} we discuss the stochastic approximation methods we use in the proofs, while the proofs themselves are in Section \ref{symmetric} for Theorem \ref{sym3main} and Section \ref{cyclic} for Theorem \ref{cyclicmain}.  In the final Section \ref{ex_sim}, we illustrate the results with some examples and simulations, including some examples beyond those covered by Theorems \ref{sym3main} and \ref{cyclicmain}.

\section{Stochastic approximation approach} \label{stochapprox}

In this section we introduce some of the stochastic approximation ideas which appear in our proofs.

For a general matrix $A$ and a given configuration of balls $x(n)$ at time $n$, let $i_{n+1} \in \{1,\ldots,d\}$ be the random colour of the ball to be added in the urn at time $n+1$.
Then, note that
\begin{align}
  \bar{x}(n+1) & = \frac{x(n) + {\bf e}_{i_{n+1}}}{n_0 + n+1} = \frac{(n_0+n )\bar{x}(n) + {\bf e}_{i_{n+1}}}{n_0+n+1} \nonumber \\
  & = \left( 1 -\frac{1}{n_0+n+1} \right)\bar{x}(n) + \frac{{\bf e}_{i_{n+1}}}{n_0+n+1},
\end{align}
implying
\begin{equation}
  \label{eq:SAA}
\bar{x}(n+1) - \bar{x}(n) = \frac{1}{n_0+n+1}({\bf e}_{i_{n+1}} - \bar{x}(n)).
\end{equation}

The idea here is to rearrange the right hand side of \eqref{eq:SAA}
into a deterministic part and a zero mean ``noise''.
More specifically, let
\begin{equation}
  \label{eq:VFF}
  F(\bar{x}(n)) := \E[ {\bf e}_{i_{n+1}} \, | \, \salj_n] - \bar{x}(n),
\end{equation}
and
\begin{equation}
  \label{eq:NOISE}
  \xi_{n+1} := {\bf e}_{i_{n+1}} - \E[ {\bf e}_{i_{n+1}} \, | \, \salj_n].
\end{equation}
By setting $\gamma_n = 1/(n_0+n+1)$, we obtain
\begin{equation}
  \label{eq:SAAX}
  \bar{x}(n+1) - \bar{x}(n) = \gamma_n(F(\bar{x}(n)) + \xi_{n+1}).
\end{equation}
The above sequence can be thought as a numerical approximating method with varying step size $\gamma_n$ for solving the ODE $dx/dt = F(x)$.
For small enough $\gamma_n$ and under mild conditions, the asymptotic behavior of $(\bar{x}_n)_{n \in \N}$ and the underlying ODE are closely connected.
This is called the \emph{ODE method} or \emph{the dynamical system approach}, which alongside some probabilistic techniques, is applied to examine almost sure dynamics of stochastic approximation processes.
In case $F$ is a gradient-like vector field and in the presence of a strict Lyapunov function whose set of critical values has empty interior,
Bena\"im \cite{Ben96, benaim} shows that the limit set $l(\bar{x}, \omega) := \bigcap_{t \geq 0}\overline{\{\bar{x}(s, \omega) \> : \> s \geq t\}}$
is almost surely a connected subset of the equilibria for the flow induced by the underlying vector field $F:\R^d\to\R^d$.

In our model, we have
$$F_i(x)=\frac{u_i}{\sum_{j=1}^N u_j}-x_i, \quad i=1,\ldots,d,$$
where $u_i$ has the same relationship to $x$ as $u_i(n)$ to $x(n)$.

\section{Proofs for the symmetric case}\label{symmetric}

In this section we prove Theorem \ref{sym3main}.

Throughout this section we let $A$ be the matrix
\begin{equation}
  \label{eq:matrix}
  \begin{pmatrix} 1 & a & a \\ a & 1 & a \\ a & a & 1 \end{pmatrix}, \quad a > 0.
\end{equation}

In this case the vector field $F$ is given by
\begin{eqnarray}F_1(x_1,x_2,x_3) &=& \frac{x_1^{\beta}+ax_2^{\beta}+ax_3^{\beta}}{(1+2a)(x_1^{\beta}+x_2^{\beta}+x_3^{\beta})}-x_1 \label{eq:1}\\
F_2(x_1,x_2,x_3) &=& \frac{ax_1^{\beta}+x_2^{\beta}+ax_3^{\beta}}{(1+2a)(x_1^{\beta}+x_2^{\beta}+x_3^{\beta})}-x_2 \label{eq:2}\\
F_3(x_1,x_2,x_3) &=& \frac{ax_1^{\beta}+ax_2^{\beta}+x_3^{\beta}}{(1+2a)(x_1^{\beta}+x_2^{\beta}+x_3^{\beta})}-x_3 \label{eq:3}\end{eqnarray}

We will first prove almost sure convergence, and then characterise the possible limits.

\begin{lem}\label{point}
  The limit set of the process $(\bar{x}(n))_{n\in\N}$, defined in \eqref{eq:proportion} with given matrix \eqref{eq:matrix}, will, almost surely, be a single point which is a stationary point of $F$.
\end{lem}
\begin{proof}
  Let
  \begin{equation}
    \label{eq:lyapunov}
    L(x_1,x_2,x_3) = -(x_1+x_2+x_3) + \frac{1}{2a+1} \left[ a\log(x_1x_2x_3) - \frac{1}{\beta}(a-1)\log(x_1^{\beta}+x_2^{\beta}+x_3^{\beta}) \right].
  \end{equation}

  Then $L$ is a strict Lyapunov function for $F$.
  In fact, straightforward differentiation gives $$\frac{\partial L}{\partial x_i} = \frac{1}{x_i} F_i.$$
  Then, denoting an integral curve of $F$ by $x(t) = (x_1(t), x_2(t), x_3(t))$, we obtain
  \[ \frac{d (L \circ x)}{d t} = \sum_{i=1}^3\frac{\partial L }{\partial x_i}\frac{d x_i}{d t} = \sum_{i=1}^3x_i\left( \frac{\partial L}{\partial x_i}\right)^2 \geq 0,\]
where the equality holds in the above inequality if and only if $F(x)=0$.

  From standard results on stochastic approximation in the presence of a strict Lyapunov function, for example Corollary 6.6 of Bena\"{i}m \cite{benaim}, the limit set of $(\bar{x}(n))_{n\in\N}$ will almost
  surely be a connected set of stationary points of $F$.
 In this case, $F$ has no connected sets of stationary points other than single points, so the limit set must be a single point, and these points are the stationary points of $F$.
\end{proof}

Note that this Lyapunov function generalises in an obvious way to more than three types, as long as all off-diagonal entries are equal.

We now need to investigate the stationary points of the vector field $F$.  Noting that the lines $x_1=x_2$, $x_1=x_3$ and $x_2=x_3$ are each invariant under $F$, define the function
\begin{equation}
  \label{eq:pol}
\mathcal{P}_{a,\beta}(z) = az^{\beta+1} - z^{\beta} + (1+a)z -2a,
\end{equation}
which we will see is related to the dynamics restricted to one of these lines. The following result shows that all stationary points of $F$ are located on at least one of these lines and expresses them in terms of solutions to $\mathcal{P}_{a,\beta}(z)=0$.

\begin{prop}\label{eqcoord}
All stationary points $(x_1,x_2,x_3)$ of $F$ have at least two of $x_1, x_2, x_3$ equal, and are of one of the forms $(\frac{r}{r+2},\frac{1}{r+2},\frac{1}{r+2})$, $(\frac{1}{r+2},\frac{r}{r+2},\frac{1}{r+2})$ or $(\frac{1}{r+2},\frac{1}{r+2},\frac{r}{r+2})$, with  $r$ a solution of $\mathcal{P}_{a,\beta}(z)=0$ in $\R^+$.

Furthermore, there are at most three possible values of $r$, one of which is always $1$, corresponding to the stationary point $(\frac13,\frac13,\frac13)$.
\end{prop}

\begin{proof}
We start off by showing that any stationary point has at least two co-ordinates equal.  We do this by writing the stationary point in the form $(x,rx,sx)$ and showing that one of $r=1$, $s=1$ or $r=s$ must hold.

Rearranging \eqref{eq:1}, \eqref{eq:2} and \eqref{eq:3} at $(x,rx,sx)$ gives \begin{eqnarray}\label{rex} x &=& \frac{1+a(r^{\beta}+s^{\beta})}{(2a+1)(1+r^{\beta}+s^{\beta})} \\  \label{rey} rx &=& \frac{r^{\beta}+a(1+s^{\beta})}{(2a+1)(1+r^{\beta}+s^{\beta})} \\ \label{rez} sx &=& \frac{s^{\beta}+a(1+s^{\beta})}{(2a+1)(1+r^{\beta}+s^{\beta})}.\end{eqnarray}
It follows that \begin{eqnarray}\label{req} r^{\beta}+a(1+s^{\beta})=r(1+a(r^{\beta}+s^{\beta})) \\ \label{seq} s^{\beta}+a(1+r^{\beta})=s(1+a(r^{\beta}+s^{\beta})).\end{eqnarray}
Take the linear combination $(s+\frac{1}{a})\times$\eqref{req}$-(r+1)\times$\eqref{seq}.  This eliminates $s^{\beta}$ and $s^{\beta+1}$, giving \begin{equation}\label{rearr3}(r^{\beta}+a)s-ar(1+r^{\beta})=(ar^{\beta}+1)s-\left(a(1+r^{\beta})+\frac{r}{a}-\frac{r^{\beta}}{a}+r^{\beta+1}-1\right),\end{equation} which can be rearranged to give \begin{equation}\label{rearr4}s(r^{\beta}-1)(a-1)=(r^{\beta+1}-1)(1-a)+(r^{\beta}-r)\left(a-\frac{1}{a}\right).\end{equation}  Assuming $a\neq 1$, \eqref{rearr4} gives $r=1$ or
\begin{equation}\label{sfromr}s=\frac{a(r^{\beta}+1)(1-r)+r^{\beta}-r}{a(r^{\beta}-1)}=\frac{\mathcal{P}_{a,\beta}(r)-ar^{\beta}+a}{a(1-r^{\beta})}.\end{equation}

Using this form for $s$ in \eqref{seq} gives (if $r\neq 1$) \begin{equation}\label{sbfromr}s^{\beta}=\frac{-ar^{\beta+1}+r^{\beta}+a-r}{a(r-1)}=\frac{\mathcal{P}_{a,\beta}(r)}{a(1-r)}+1.\end{equation}  Combining \eqref{sfromr} and \eqref{sbfromr} tells us that either $r=1$ or $s=1$ or $$\frac{s^{\beta}-1}{s-1}=\frac{r^{\beta}-1}{r-1},$$ and the latter case implies $r=s$.  Hence any stationary point has two co-ordinates equal.

We now assume, without loss of generality, that the stationary point is of the form $(x,x,rx)$ or equivalently $\left(\frac{1}{r+2},\frac{1}{r+2},\frac{r}{r+2}\right)$.  That the stationary point equations for a point of this form imply $\mathcal{P}_{a,\beta}(r)=0$ is easy to check, and it is also easy to check that $\mathcal{P}_{a,\beta}(1)=0$ for any $a,\beta>0$.

The function $\mathcal{P}_{a,\beta}$ satisfies $\mathcal{P}_{a,\beta}(0)<0$ and $\mathcal{P}_{a,\beta}(z)\to \infty$ as $z\to\infty$; furthermore it is concave for $z<\frac{\beta-1}{a(\beta+1)}$ and convex for $z>\frac{\beta-1}{a(\beta+1)}$, which indicates that it has either one root or three in $\R^{+}$, counting multiplicity.  This completes the proof.

\end{proof}

\begin{prop}\label{Proots} Consider $\mathcal{P}_{a,\beta}(z) = az^{\beta+1} - z^{\beta} + (1+a)z -2a$ for $z\in \R^+$, with $a,\beta > 0$.
  \begin{enumerate}
  \item For a given value of $a>1$, $\mathcal{P}_{a,\beta}$ has only one root at $z=1$.\label{a>1}
  \item For a given value of $a<1$, there exists $\beta_1(a)$ satisfying $\frac{1+2a}{1-a}>\beta_1(a)>1$ such that
    \begin{enumerate}
      \item If $\beta>\frac{1+2a}{1-a}$ then $\mathcal{P}_{a,\beta}'(1)<0$ and we have that $\mathcal{P}_{a,\beta}$ has three roots in $\R^+$, $1$, $r_{-}$ and $r_{+}$, labelled so that $r_{-}<1<r_{+}$.  As functions of $\beta$ for fixed $a$, $r_{+}$ is increasing and $r_{-}$ is decreasing.
      \item If $\beta_1(a)<\beta<\frac{1+2a}{1-a}$ then $\mathcal{P}_{a,\beta}'(1)>0$ and $\mathcal{P}_{a,\beta}$ has three roots in $\R^+$, $1$, $r_1$ and $r_2$, labelled so that $1<r_1<r_2$.  As functions of $\beta$ for fixed $a$, $r_{2}$ is increasing and $r_{1}$ is decreasing.
      \item If $\beta<\beta_1(a)$ then $\mathcal{P}_{a,\beta}'(1)>0$ and the only root of $\mathcal{P}_{a,\beta}$ in $\R^+$ is $1$.
    \end{enumerate}
    Furthermore $\beta_1(a)$ is an increasing function of $a$ with $\beta_1(a) \to\infty$ as $a\to 1$.
  \end{enumerate}
\end{prop}
\begin{proof}
  \emph{(i)} First, if $\beta<1$ we have that $\mathcal{P}_{a,\beta}''(z)>0$ in $\R^+$ implying that $\mathcal{P}_{a,\beta}'(z)$ is strictly increasing.
  Moreover, $\mathcal{P}_{a,\beta}'(z)$  goes from $-\infty$ to $+\infty$ when $z$ ranges from $0$ to $+\infty$ and $\mathcal{P}_{a,\beta}'(1) > 0 $.
  Then $\mathcal{P}_{a,\beta}'(z)$ changes sign only once at some $z^*<1$.
  Now, since $\mathcal{P}_{a,\beta}(0) = -2a < 0$ and $\mathcal{P}_{a,\beta}(z)$ is decreasing for $z<z^*<1$ and increasing otherwise, it follows that $\mathcal{P}_{a,\beta}(z)$ crosses the $z=0$ line only once at $z=1$.
  Second, the same happens for $\beta > 1$ since $\mathcal{P}_{a,\beta}'(z) > 0$ in $\R^+$ and so $\mathcal{P}_{a,\beta}(z)$ is strictly increasing.
  The case $\beta = 1 $ is trivial.

  \emph{(ii)} If $\beta>\frac{1+2a}{1-a}$ we have $\mathcal{P}_{a,\beta}'(1)<0$, indicating that in this case $\mathcal{P}_{a,\beta}$ must have three roots.
  Note that if $\beta=\frac{1+2a}{1-a}$ then $\mathcal{P}_{a,\beta}'(1)=0$ but that $\mathcal{P}_{a,\beta}''(1)<0$, showing that this is a double root, not a triple root, and so there must be another root in that case for larger $z$.
  If $\beta<\frac{1+2a}{1-a}$ then  $\mathcal{P}_{a,\beta}'(1)>0$ and $\mathcal{P}_{a,\beta}$ has no root less than $1$.
  Then, there must be either none, one double, or two distinct additional roots greater than 1.

The derivative of $\mathcal{P}_{a,\beta}(z)$ with respect to $\beta$, for fixed $a$ and $z$, is $(az-1)(\log z)z^{\beta}$.  Thus, for any $z\in(1,1/a)$, $\mathcal{P}_{a,\beta}(z)$ is decreasing in $z$, and it follows that if there are roots of $\mathcal{P}_{a,\beta}$ in this range for a particular value of $\beta$ there must also be for any larger $\beta$.  As $\mathcal{P}_{a,\beta}(z)>0$ if $z\geq 1/a$, this shows that there exists $\beta_1(a)\in[1,\frac{1+2a}{1-a}]$ such that there is one root of $\mathcal{P}_{a,\beta}$ when $\beta<\beta_1(a)$ and three when $\beta>\beta_1(a)$.

Let $\beta_0(a)>\frac{2}{1-a}-1>1$ be the unique solution to \begin{equation}\label{beta0}1+a=\left(\left(1-\frac{2}{\beta+1}\right)\frac{1}{a}\right)^{\beta-1}.\end{equation} (It can be seen that \eqref{beta0} has a unique solution for fixed $a<1$, as in that case the right hand side is increasing in $\beta$ if the right hand side is greater than $1$, is equal to $1$ at $\beta=1$ and tends to $\infty$ as $\beta\to\infty$.  That $\beta_0(a)>\frac{2}{1-a}-1$ can be seen by noting that if $\beta$ is a solution of \eqref{beta0} we must have $\left(1-\frac{2}{\beta+1}\right)\frac{1}{a}>1$.)   Then if $\beta<\beta_0(a)$ we have $\mathcal{P}_{a,\beta}'(z)>0$ for all $z$ and hence $\mathcal{P}_{a,\beta}$ is increasing in $z$ and so $z=1$ is the only root.  This shows that $\beta_1(a) \geq \beta_0(a)>1$.

Now, the fact that as mentioned above there is a root greater than $1$ when $\beta=\frac{1+2a}{1-a}$ together with the continuity of $\mathcal{P}_{a,\beta}(z)$ in $\beta$ ensures that there remains a root greater than $1$ for $\beta\in\left(\frac{1+2a}{1-a}-\epsilon,\frac{1+2a}{1-a}\right)$ for some $\epsilon>0$, so $\beta_1(a)<\frac{1+2a}{1-a}$.

The claims that $r_{+}$ and $r_2$ are increasing functions of $\beta$ and that $r_1$ is a decreasing function of $\beta$ also follow from the negative derivative of $\mathcal{P}_{a,\beta}(z)$ with respect to $\beta$ for $z\in(1,1/a)$.  Similarly the claim that $r_{-}$ is a decreasing function of $\beta$ follows from the derivative of $\mathcal{P}_{a,\beta}(z)$ with respect to $\beta$ being positive on $(0,1)$.

To see that  $\beta_1(a)$ is an increasing function of $a$, note that for fixed $z>1$ the derivative of $\mathcal{P}_{a,\beta}(z)$ with respect to $a$ is positive, meaning that if we are in case (c) for particular choices of $a$ and $\beta$ we will also be in case (c) for the same value of $\beta$ and any larger value of $a$.  That $\beta_1(a)\to \infty$ as $a\to 1$ follows from $\beta_0(a)>\frac{2}{1-a}-1$.
\end{proof}

We now investigate the stability of these roots. We will consider a stationary point $p$ of $F$ to be stable if it is a local maximum of the Lyapunov function $L$, and to be unstable if it is a saddle point or local minimum of $L$.  Furthermore, if all eigenvalues of $DF(p)$ have negative real part, $p$ is said to be linearly stable, while if some eigenvalue has positive real part, $p$ is said to be linearly unstable.

\begin{prop}\label{stable}
If a stationary point for $F$ is of the form $(x,x,rx)$ or $(x,rx,x)$ or $(rx,x,x)$, then it is linearly stable if $\mathcal{P}_{a,\beta}'(r)>0$ and $\frac{r^{\beta}+2}{r+2}>\frac{\beta(1-a)}{2a+1}$, and linearly unstable if either $\mathcal{P}_{a,\beta}'(r)<0$ or $\frac{r^{\beta}+2}{r+2}<\frac{\beta(1-a)}{2a+1}$.\end{prop}

\begin{proof}
Without loss of generality we focus on the case $(x,x,rx)$.

We note that the differential equation driven by $F$ keeps the line $x_1=x_2$ invariant, so we consider it restricted to this line;  the equation for $F_3$ gives $$F_3\left(\frac{1-x_3}{2},\frac{1-x_3}{2},x_3\right)=\frac{2a\left(\frac{1-x_3}{2}\right)^{\beta}+x_3^{\beta}}{(1+2a)\left(2\left(\frac{1-x_3}{2}\right)^{\beta}+x_3^{\beta}\right)}-x_3.$$  Let $x_3=z/(z+2)$ so that $x_1=x_2=1/(z+2)$.  Then $$F_3\left(\frac{1-x_3}{2},\frac{1-x_3}{2},x_3\right)=\frac{-2\mathcal{P}_{a,\beta}(z)}{(1+2a)(2+z^{\beta})},$$ and so is positive when $\mathcal{P}_{a,\beta}(z)$ is negative and vice versa.  Hence a stationary point $(x,x,rx)$ is stable in this direction if $\mathcal{P}_{a,\beta}'(r)>0$ and unstable in this direction if $\mathcal{P}_{a,\beta}'(r)<0$.

Because $F$ is symmetric in $x_1$ and $x_2$, the other direction in which we need to consider stability will be perpendicular to this one.  Hence we consider
\begin{equation}\label{estability}
  \begin{split}
  F_1(x+\epsilon,x-\epsilon,rx) &= \frac{(x+\epsilon)^{\beta}+(x-\epsilon)^{\beta}+ar^{\beta}x^{\beta}}{(2a+1)((x+\epsilon)^{\beta}+(x-\epsilon)^{\beta}+r^{\beta}x^{\beta})}-x-\epsilon \\
                                &= F_1(x,x,rx)+\epsilon\left(-1+\frac{\beta(1-a)}{(2a+1)(2+r^{\beta})x}\right)+o(\epsilon) \\
                                &= F_1(x,x,rx)+\epsilon\left(-1+\frac{\beta(1-a)(r+2)}{(2a+1)(2+r^{\beta})}\right)+o(\epsilon).
  \end{split}
\end{equation}

It follows that $(x,x,rx)$ is a stable stationary point in the direction perpendicular to the line $x_1=x_2$ if $\frac{r^{\beta}+2}{r+2}>\frac{\beta(1-a)}{2a+1}$ and unstable in that direction if the reverse inequality applies.
\end{proof}

We shall henceforth restrict ourselves to the case $a<1$ since Propositions \ref{eqcoord}, \ref{Proots}(\ref{a>1}) and \ref{stable} imply that if $a>1$, $(\frac13,\frac13,\frac13)$ is the only stable stationary point for $F$ and by Corollary \ref{point}, $(\bar{x}(n))_{n \in \N}$ must converge to it.  The case $a=1$ has the probabilities of each colour being $1/3$ regardless of $\bar{x}(n)$ and so it is easily seen that $\bar{x}(n)\to (1/3,1/3,1/3)$ almost surely.

\begin{cor}\label{criteria}
  Assume $a<1$.
  \begin{enumerate}

  \item If $\beta<\beta_1(a)$, then the stationary point $(\frac13,\frac13,\frac13)$ is stable, and is the limit with probability $1$.

  \item If $\beta_1(a)<\beta<\frac{1+2a}{1-a}$, then the stationary points $(\frac13,\frac13,\frac13)$ and $\left(\frac{1}{r_{2}+2},\frac{1}{r_{2}+2},\frac{r_2}{r_{2}+2}\right)$ (and its permutations) are  stable, while the stationary point $\left(\frac{1}{r_{1}+2},\frac{1}{r_{1}+2},\frac{r_1}{r_{1}+2}\right)$ and its permutations are linearly unstable.

  \item If $\beta>\frac{1+2a}{1-a}$, then there are three stationary points of $F$ of the form $(x,x,rx)$ corresponding to the three solutions $r_{-}<1<r_{+}$ of $\mathcal{P}_{a,\beta}(z)=0$ in $\R^+$.  The stationary points $(\frac13,\frac13,\frac13)$ and $\left(\frac{1}{r_{-}+2},\frac{1}{r_{-}+2},\frac{2}{r_{-}+2}\right)$ (and its permutations) are linearly unstable, while $\left(\frac{1}{r_{+}+2},\frac{1}{r_{+}+2},\frac{r_{+}}{r_{+}+2}\right)$ and its permutations are stable.

  \end{enumerate}

\end{cor}

\begin{proof}

  \emph{(i)}
  Stability follows from Proposition \ref{stable}, and almost sure convergence from Corollary \ref{point}.
  \emph{(ii)}
  The shape of $\mathcal{P}_{a,\beta}$ ensures that $r_1,r_2>1$ and that $\mathcal{P}_{a,\beta}'(r_1)<0$ and $\mathcal{P}_{a,\beta}'(r_2)>0$, showing that $\left(\frac{1}{r_{1}+2},\frac{1}{r_{1}+2},\frac{r_1}{r_{1}+2}\right)$ is unstable, and that for the other two stationary points we just need to check the stability perpendicular to the line $x_1=x_2$.  But $\frac{r_2^{\beta}+2}{r_2+2}>1>\frac{\beta(1-a)}{2a+1}$ by our assumption on $\beta$, so the condition from Propostion \ref{stable} is satisfied and so $(\frac13,\frac13,\frac13)$ and $\left(\frac{1}{r_{2}+2},\frac{1}{r_{2}+2},\frac{r_2}{r_{2}+2}\right)$ are stable.
    \emph{(iii)}
  That $\mathcal{P}_{a,\beta}(z)$ has three solutions follows from Proposition \ref{Proots}. As $\mathcal{P}_{a,\beta}'(1)<0$ it follows from Propostion \ref{stable} that $(\frac13,\frac13,\frac13)$ is unstable, and as $r_{-}<1$ we have $\frac{r_{-}^{\beta}+2}{r_{-}+2}<1<\frac{\beta(1-a)}{2a+1}$, so $\left(\frac{1}{r_{-}+2},\frac{1}{r_{-}+2},\frac{r_{-}}{r_{-}+2}\right)$ is also unstable.  The global maximum of the Lyapunov function on $\Delta^2$ must be a stable stationary point of $F$, so the remaining stationary points, $\left(\frac{1}{r_{+}+2},\frac{1}{r_{+}+2},\frac{r_{+}}{r_{+}+2}\right)$ and its permutations, must be stable.
\end{proof}

Since our process $(\bar{x}(n))_{n \in \N}$ is in some sense a discrete version of the differential flow $dx(t)/dt = F(x(t))$, we would like to ascertain whether the stochastic process may or may not converge to certain critical points in terms of what type of critical points they are for the associated flow.

A subset $A$ of the phase state is called an attractor for a flux $\Phi$ if it is a nonempty, compact and invariant subset having a  neighbourhood $W$ such that the distance $d(\Phi_t, A) \rightarrow 0$ as $t \rightarrow \infty$ uniformly in $W$.
Now, let $p$ be a stationary point of a smooth vector field $F$.  If $p$ is a stable fixed point of $F$, then $p$ is an attractor for the flux induced by $F$.

The following result completes the proof of Theorem \ref{sym3main}.

\begin{prop}\label{convergence} Let $l(\bar{x})$ denote the limit set of the process $(\bar{x}(n))_{n \in \N}$ defined in \eqref{eq:SAAX} and recall the stability criteria in Proposition \ref{stable}. Then we have
  \begin{enumerate}
  \item $\pr[l(\bar{x}) = \{p\}] > 0$ for stable points $p$ of $F$.
  \item $\pr[l(\bar{x}) = \{p\}] = 0$ for linearly unstable points $p$ of $F$.
  \end{enumerate}
\end{prop}

\begin{proof}
  Without loss of generality we focus on the case $(x,x,rx)$.

  \emph{(i)} Let us now show that the process $\bar{x}$ in fact converges with positive probability toward a given attractor.
  Of course, it is necessary that the process has positive probability of being arbitrarily close to the attractor at arbitrarily large times.
  That is, a point $p$ is said to be attainable by a process $X$ if for each $t>0$ we have that  $\pr[\exists ~ s \geq t ~:~ X(s) \in N_p] > 0$ for every neighborhood $N_p$ of $p$.
  It turns out that if the function $F + Id$ associated with an urn process $X$ maps the simplex into its interior, it follows that every point of the simplex is attainable by $X$.
  This is indeed the case for our process $\bar{x}$.
  Finally, Bena\"{i}m \cite{benaim} Theorem 7.3 ensures that a given attainable attractor $p$ with non-empty basin of attraction is such that $\pr[l(\bar{x}) = \{p\}] > 0$.

  \emph{(ii)} Let $p$ be a linearly unstable critical point and $\mathcal{N}_p \subset \Delta_2$ a neighborhood of $p$.
  The simplex is considered as a differential manifold by identifying its tangent space at any point with the linear subspace
  $T\Delta_2 = \{x \in \R^3 \> : \> \sum_i x_i = 0\}$.
  In our case, the only non-trivial condition ensuring Pemantle's non-convergence criteria (see \cite{pemantle1990}) is that we have positive expectation of the positive part of the component of the noise in any given direction.
  Formally, we need that whenever $\bar{x}_n \in \mathcal{N}_p$, there is a constant $\kappa$ such that
  $\E[\max\{\xi_{n+1} \cdot \theta,0\}  \> | \> \mathcal{F}_n] \geq \kappa$
  for every unit vector $\theta \in T\Delta_2$.
  For notational simplicity, write $\tilde{u}_i = u_i(n)/(\sum_j u_j(n))$ and note that
\begin{align}\label{eq:drift}
  \E[\max\{\xi_{n+1} \cdot \theta,0\}  \> | \> \mathcal{F}_n]  =
  & \tilde{u}_1 \max\{\theta_1(1-\tilde{u}_1)-\theta_2 \tilde{u}_2-\theta_3 \tilde{u}_3,0\}  + \nonumber \\
  & \tilde{u}_2 \max\{-\theta_1 \tilde{u}_1+\theta_2(1-\tilde{u}_2)-\theta_3 \tilde{u}_3),0\} + \nonumber \\
  & \tilde{u}_3 \max\{\theta_1 \tilde{u}_1-\theta_2 \tilde{u}_2+\theta_3(1-\tilde{u}_3),0\}.
\end{align}
Now, write $\theta = (\theta_1,\theta_2,\theta_3)$ and suppose that $\theta_i < 0$ for exactly two coordinates.
Then, it is simple to find a positive term in \eqref{eq:drift}.
On the other hand, suppose that $\theta_i < 0$ for exactly one coordinate.
Let us say $\theta_3 < 0$.
Then, depending on which of the inequalities $\theta_1 > \theta_2 \geq  0$ or $\theta_2 > \theta_1 \geq  0$ holds, it is also not difficult to find a positive term in \eqref{eq:drift}.
Finally, to prove that it is in fact uniformly positive, it is enough that $\theta$ is an unit vector in $T\Delta_2$ and that all stationary points $p$ are such that $\tilde{u}_i$ is uniformly positive in a neighborhood $\mathcal{N}_p$ of $p$ for $i=1,2,3$.
\end{proof}

\section{Proofs for the cyclic case}\label{cyclic}

\subsection{Introduction}

In this section we prove Theorem \ref{cyclicmain}. Here $A$ is the matrix
$$\begin{pmatrix} 1 & 1 & 0 \\ 0 & 1 & 1 \\ 1 & 0 & 1 \end{pmatrix},$$ and we have \begin{eqnarray*}F_1(x_1,x_2,x_3) &=& \frac{x_1^{\beta}+x_2^{\beta}}{2(x_1^{\beta}+x_2^{\beta}+x_3^{\beta})}-x_1\\
F_2(x_1,x_2,x_3) &=& \frac{x_2^{\beta}+x_3^{\beta}}{2)(x_1^{\beta}+x_2^{\beta}+x_3^{\beta})}-x_2\\
F_3(x_1,x_2,x_3) &=& \frac{x_3^{\beta}+x_1^{\beta}}{2(x_1^{\beta}+x_2^{\beta}+x_3^{\beta})}-x_3\end{eqnarray*}

First, we note that for any choice of $\beta$, $(\frac13,\frac13,\frac13)$ is a stationary point of $F$.
The following two results give information on its stability and show that it is in fact the only stationary point.

\begin{lem}\label{stability}For the vector field $F$, the stationary point at $(\frac13,\frac13,\frac13)$ is stable if $\beta<4$, and a linearly unstable source if $\beta>4$. \end{lem}

\begin{proof}
As we are working with $x\in\Delta^2$, write $x_3=1-x_1-x_2$.  Routine calculus then shows that the Jacobian matrix at $(\frac13,\frac13,\frac13)$ is $$\begin{pmatrix}\frac{\beta}{2}-1 & \frac{\beta}{2} \\ -\frac{\beta}{2} &  -1 \end{pmatrix}.$$
The eigenvalues of this Jacobian are then the roots $\lambda$ of $$\lambda^2+\lambda\left(2-\frac{\beta}2\right)-\left(\frac{\beta}2-1\right)+\left(\frac{\beta}{2}\right)^2$$
which are $$\left(\frac{\beta}{4}-1\right)\pm \frac{i\sqrt{3}\beta }{4}.$$  The result follows.
\end{proof}

\begin{lem}\label{statpoint}The only stationary point of $F$ in $\Delta^2$ is $(\frac13,\frac13,\frac13)$.\end{lem}

\begin{proof}
For $x=(x_1,x_2,x_3)\in \Delta^2$ with $F_1(x)=F_2(x)=F_3(x)=0$, we have $x_1=\frac{x_1^{\beta}+x_2^{\beta}}{x_2^{\beta}+x_3^{\beta}}x_2$ and similarly $x_2=\frac{x_2^{\beta}+x_3^{\beta}}{x_3^{\beta}+x_1^{\beta}}x_3$ and $x_3=\frac{x_3^{\beta}+x_1^{\beta}}{x_1^{\beta}+x_2^{\beta}}x_1$.  Using this, $$x_1-x_2=x_2\frac{x_1^{\beta}-x_3^{\beta}}{x_2^{\beta}+x_3^{\beta}},$$ indicating that (if $x_2>0$) if $x_1>x_2$ then also $x_1>x_3$, while if $x_1<x_2$ then $x_1<x_3$.  Similarly, if $x_3>0$ then the signs of $x_2-x_3$ and $x_2-x_1$ are the same, and if $x_1>0$ then the signs of $x_3-x_1$ and $x_3-x_2$ are the same.  Hence the only stationary point of $F$ in the interior of $\Delta^2$ is $(\frac13,\frac13,\frac13)$.

It is also easy to check that if $x_1=0$ then $x_2=0$, and similarly that if $x_2=0$ then $x_3=0$ and if $x_3=0$ then $x_1=0$.  Hence there are no stationary points of $F$ on the boundary of $\Delta^2$.
\end{proof}

We can now complete the proof of Theorem \ref{cyclicmain}.

That we only have one stationary point, and that it is never a saddle, restricts the possibilities for chain transitive sets.  In two dimensions Theorem 6.12 of Bena\"{i}m \cite{benaim} states that chain transitive sets must be unions of stationary points, periodic orbits and cyclic orbit chains.  However, with only one stationary point which is not a saddle cyclic orbit chains are impossible.  We can thus conclude that the limit set must be a connected union of periodic orbits and stationary points.

By Lemma \ref{stability}, if $\beta<4$ then by Lemma \ref{stability} the stationary point $(\frac13,\frac13,\frac13)$  is stable, and hence is an attractor for the flow given by $F$.  By Theorem 7.3 of Bena\"{i}m \cite{benaim}, to show that there is positive probability of convergence to $(\frac13,\frac13,\frac13)$ it is enough to show that it is an attainable point, that is that the process has positive probability of being arbitrarily close to $(\frac13,\frac13,\frac13)$ at arbitrarily large times.  This is straightforward to show: for any $\epsilon>0$ there will be points of the form $\left(\frac{n_1}{n},\frac{n_2}{n},\frac{n_3}{n}\right)$ with $n_1,n_2,n_3$ integers satisfying $n_1+n_2+n_3=n$ and $\max\left\{\frac{n_1}{n}-\frac13,\frac{n_2}{n}-\frac13,\frac{n_3}{n}-\frac13\right\}<\epsilon$ for arbitrarily large $n$.  There will be positive probability of any such point being reached, so $(\frac13,\frac13,\frac13)$ is indeed attainable.

If $\beta>4$, then by Lemma \ref{stability} $(\frac13,\frac13,\frac13)$ is linearly unstable, so as in Proposition \ref{convergence} it will follow that it is a limit with probability zero if we have positive expectation of the positive part of the component of the noise in any given direction.  In fact, the same argument as in Proposition \ref{convergence} will work here, so we can conclude that $\bar{x}(n)$ has probability zero of converging to a stationary point.  It follows that the limit set must be a periodic orbit or a connected union of periodic orbits, completing the proof of Theorem \ref{cyclicmain}.

\section{Examples and simulations} \label{ex_sim}

In this section we consider some examples, including some where exact calculations are possible, and some simulations.  We also consider some examples which go beyond the cases covered by Theorems \ref{sym3main} and \ref{cyclicmain}.

\subsection{The symmetric case with $\beta=2$}
In the case of Theorem \ref{sym3main} with $\beta=2$, the possible limits and the phase transitions can be explicitly identified.  We find that $\mathcal{P}_{a,2}(z)=(z-1)(az^2+(a-1)+2a)$, with roots given by $z=1$ and $z=\frac{1-a \pm\sqrt{1-2a-7a^2}}{2a}$.  If $a< \frac{\sqrt{8}-1}{7}$, then these are real and positive, and letting $\lambda_1=\frac{3a+1-\sqrt{1-2a-7a^2}}{4(2a+1)}, \lambda_2=\frac{3a+1+\sqrt{1-2a-7a^2}}{4(2a+1)}, \lambda_3=\frac{a+1-\sqrt{1-2a-7a^2}}{2(2a+1)}$ and $\lambda_4=\frac{a+1+\sqrt{1-2a-7a^2}}{2(2a+1)}$, we obtain linearly stable stationary points $(\lambda_1,\lambda_1,\lambda_4)$, $(\lambda_1,\lambda_4,\lambda_1)$ and $(\lambda_4,\lambda_1,\lambda_1)$, and linearly unstable stationary (except at $a=\frac14$) points $(\lambda_2,\lambda_2,\lambda_3)$, $(\lambda_2,\lambda_3,\lambda_2)$, $(\lambda_3,\lambda_2,\lambda_2)$.  In addition, the stationary point $(\frac13,\frac13,\frac13)$ is linearly stable if $a>\frac14$, and linearly unstable if $a<\frac14$.

If $a>\frac{\sqrt{8}-1}{7}$, then $(\frac13,\frac13,\frac13)$ is the only stationary point, and is stable.  In the notation of Theorem \ref{sym3main}, we have $\beta_1\left(\frac{\sqrt{8}-1}{7}\right)=2$, and the three phases are as follows:

\begin{figure}[h]\label{fig:simulations-sym}
  \caption{20 simulations of the symmetric model for $\beta=2$}
  \centering
\begin{subfigure}{.33\textwidth}
  \centering

  \includegraphics[width=5cm, height=5cm]{Rplot11}
  \caption{$a=0.2$}\label{ap2}
  \label{fig:sub1}
\end{subfigure}
 \centering
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=5cm, height=5cm]{Rplot22}
  \caption{$a=0.26$}\label{ap26}
  \label{fig:sub2}
\end{subfigure}
\centering
\begin{subfigure}{.33\textwidth}
  \centering

  \includegraphics[width=5cm, height=5cm]{Rplot33}
  \caption{$a=0.5$}\label{ap5}
  \label{fig:sub3}
\end{subfigure}
\end{figure}

\begin{itemize}
\item When $a<\frac14$, $(\frac13,\frac13,\frac13)$ is unstable; there are three other stationary points, $(\lambda_1,\lambda_1,\lambda_4)$ and permutations, placed symmetrically, which are stable.  For example when $a=\frac15$, $(0.1847,0.1847,0.6306)$ and permutations are stable.  Almost surely, one of these three points will be the limit.  A simulation of 20 trajectories of the stochastic process in this case appears in Figure \ref{ap2}.
\item  For $\frac14<a<\frac{\sqrt{8}-1}{7}$, $(\frac13,\frac13,\frac13)$ is now stable but there are also stable stationary points elsewhere, near $(\frac14,\frac14,\frac12)$.  In this case, both symmetric and asymmetric limits have positive probability.  For example, at $a=0.26$, there are stable stationary points at $(0.2792,0.2792,0.4416)$ and permutations as well as $(\frac13,\frac13,\frac13)$.  A simulation of 20 trajectories of the stochastic process in this case appears in Figure \ref{ap26}.
\item For $a>\frac{\sqrt{8}-1}{7}$, $(\frac13,\frac13,\frac13)$ is the only stationary point, and is stable, and will be the limit almost surely.  A simulation of 20 trajectories of the stochastic process in the $a=\frac12$ case appears in Figure \ref{ap5}.
\end{itemize}

At the critical value $a=\frac14$, $(\frac13,\frac13,\frac13)=(\lambda_2,\lambda_2,\lambda_3)$ has zeros as eigenvalues of its Jacobian and so is neither linearly stable nor linearly unstable, while there are stable stationary points at $(\frac14,\frac14,\frac12)$ and permutations; similarly at the critical value $a=\frac{\sqrt{8}-1}{7}$ the stationary point $(\lambda_2,\lambda_2,\lambda_3)=(\lambda_1,\lambda_1,\lambda_4)$ is neither linearly stable nor linearly unstable.

\subsection{The symmetric case with $\beta=3$}
It is also possible to do some explicit calculations when $\beta=3$.  In this case $(\frac13,\frac13,\frac13)$ is linearly stable when $a>\frac25$ and linearly unstable when $a<\frac25$, and we have $$\mathcal{P}_{a,3}(z)=(z-1)(az^3+(a-1)z^2+(a-1)z+2a),$$ where the cubic factor has one real root (which is negative) when $a>a_c=\frac{1}{166}(3.(2)^{1/2}.(3)^{1/4} + 24.(3)^{1/2} + 13.(2)^{1/2}.(3)^{3/4}-20) = 0.4160306$ and three real roots (one of which is negative) when $a<a_c$.  (In the notation of Theorem \ref{sym3main}, $\beta_1(a_c)=3$.)  Hence for $a>a_c$ we get almost sure convergence to $(\frac13,\frac13,\frac13)$, for $a<\frac25$ we get almost sure convergence to one of three asymmetric stationary points, and for $\frac25<a<a_c$ the process may converge either to $(\frac13,\frac13,\frac13)$ or to an asymmetric stationary point, each with positive probability.

\subsection{The cyclic model}

Figure \ref{cyclic3} shows 20 simulations of the cyclic model when $\beta=3$, showing convergence to  $(\frac13,\frac13,\frac13)$.
Figure \ref{cyclic6} shows 20 simulations with $\beta=6$, showing apparent convergence to a single limit cycle.

\begin{figure}[h]\label{fig:simulations-cyclic}
  \caption{20 simulations of the cyclic model}
  \centering
\begin{subfigure}{.33\textwidth}
  \centering

  \includegraphics[width=5cm, height=5cm]{Rplot44}
  \caption{$\beta=3$}\label{cyclic3}
  \label{fig:sub1-cyclic}
\end{subfigure}
 \centering
\begin{subfigure}{.33\textwidth}
  \centering
  \includegraphics[width=5cm, height=5cm]{Rplot55}
  \caption{$\beta=6$}\label{cyclic6}
  \label{fig:sub2-cyclic}
\end{subfigure}
\end{figure}

\subsection{The symmetric case with more than three types}

It is natural to extend Section \ref{symmetric} to $d>3$ types, letting $A$ be the $d\times d$ matrix with $a_{ii}=1$ for each $i$ and $a_{ij}=a$ for $i\neq j$.  It is straightforward to extend the Lyapunov function \eqref{eq:lyapunov} to this case, meaning that Lemma \ref{point} applies.  However, the later calculations, starting with Proposition \ref{eqcoord}, do not apply.  It thus may be interesting to investigate whether more complex patterns of phases can occur in this case than when $d=3$; however, numerical solution of the stationary point equations for particular values of $a$ when $d=4,5,6$ suggests that the behaviour is in fact very similar to the $d=3$ case, with three phases which parallel those found in Theorem \ref{sym3main}.

\subsection{A more general cyclic case}

It is natural to extend Section \ref{cyclic} by allowing the matrix $A$ to take the form $$\begin{pmatrix} 1 & a & 0 \\ 0 & 1 & a \\ a & 0 & 1\end{pmatrix},$$ allowing for different strengths of the cyclic reinforcement.  It is straightforward to extend Lemma \ref{stability} to this case, showing that the stationary point at $(\frac13,\frac13,\frac13)$ is linearly stable if $a\geq 2$ or if $a<2$ and $\beta<\frac{2(1+a)}{2-a}$, and linearly unstable if $a<2$ and $\beta>\frac{2(1+a)}{2-a}$.  However Lemma \ref{statpoint} does not apply for general $a$ and there may be other stationary points.

Numerical investigation when $\beta=2$ suggests that there are three phases in $a$: in addition to a phase with almost sure convergence to $(\frac13,\frac13,\frac13)$ and a phase with convergence to a limit cycle, there is a phase up to $a \approx 0.25057$ where there are stable stationary points other than $(\frac13,\frac13,\frac13)$ and that the process usually converges to one of these.

\section*{Acknowledgements}

M.C.’s research was supported by CONICET [grant number 10520170300561CO].
The authors are grateful to Andrew Wade for suggestions that improved the presentation.

\bibliographystyle{plain}
\bibliography{references}

\end{document}

config: arxiv


text:
\section{Conclusion}\label{sec:summary-future-directions}
We proposed a confidence sequence for bounded random variables and
examined its efficacy in synthetic examples and applications.
Our method is favorable to methods that utilize classical concentration results. It can be applied to various problems for improved performance, including testing equality of distributions, testing independence~\cite{balsubramani2016sequential}, etc.
Our work can be extended in a few future directions. We assumed that $X_i$'s are independent and bounded. The generalizations for the dependent case and/or the sub-Gaussian case are of interest. The generalized results can be obtained based on the results of \citet[Theorem 2.1]{pinelis2006binomial} and \citet{bentkus2010bounds}.

Regarding dependence, Theorem~2.1 of~\citet{pinelis2006binomial} shows that Theorem~\ref{thm:initial-maximal-ineq} holds even if $X_i$'s form a supermartingale difference sequence
, i.e., assumption~\eqref{eq:boundedness-conditions} is replaced by $\mathbb{E}[X_i|X_1, \ldots, X_{i-1}] \le 0$, $\mathbb{P}(X_i > B) = 0$, and
$\mathbb{P}(\mathbb{E}[X_i^2|X_1, \ldots, X_{i-1}] \le A_i^2) = 1.$
Theorem~\ref{thm:uniform-in-n} follows readily, but Theorem~\ref{thm:uniform-in-n-empirical-Bentkus} requires further restrictions that allow estimation of $A^2_i$.

The boundedness assumption~\eqref{eq:boundedness-conditions}, which maybe restrictive for applications in statistics, finance and economics,
can be replaced by
\begin{equation}\label{eq:generalize-bnd-conditions}
\textstyle
\mathbb{E}[X_i] = 0,\;\mbox{Var}(X_i) \le A_i^2,\quad\mbox{and}\quad \mathbb{P}(X_i > x) \le \widebar{F}(x),
\end{equation}
for all $x\in\mathbb{R}$, where $\widebar{F}(\cdot)$ is a survival function on $[0, \infty)$, i.e., $\widebar{F}(\cdot)$ is non-increasing and $\widebar{F}(0) = 1$ and $\widebar{F}(\infty) = 0$.
For example,
$\widebar{F}(x) = 1/\{1 + (x/K)_+^{\alpha}\}$, {or} $\widebar{F}(x) = \exp(-(x/K)_+^{\alpha})$, {for some} $K > 0$;
$\alpha = 2$ in the second example is sub-Gaussianity.
Akin to~\eqref{eq:main-optimality}, there exist random variables $\eta_i = \eta_i(A_i^2, \widebar{F})$ satisfying~\eqref{eq:generalize-bnd-conditions} such that
$$\textstyle \sup_{X_i\sim~\eqref{eq:generalize-bnd-conditions}}\mathbb{E}[\left(\sum_{i=1}^n X_i - t\right)_+^2] = \mathbb{E}[\left(\sum_{i=1}^n \eta_i - t\right)_+^2], $$
for all $n\ge1$ and $t\in\mathbb{R}$. Here the superemum is taken over all distributions satisfying~\eqref{eq:generalize-bnd-conditions}. Hence, Theorem~\ref{thm:initial-maximal-ineq} can be generalized, which in turn leads to generalizations of Theorems~\ref{thm:uniform-in-n} and~\ref{thm:uniform-in-n-empirical-Bentkus}. The details on the construction of~$\eta_i$ and the corresponding confidence sequence will be discussed elsewhere.


TITLE: Why can the retarded propagator be defined in terms of a commutator or without it?
QUESTION [3 upvotes]: In many textbooks [e.g. Peskin & Schröder p. 30 eq. (2.55), or Tong's notes p. 41 eq. (2.101)], the retarded propagator is defined as
$$G_R = \Theta(x^0-y^0) \left< [\phi(x), \phi(y)] \right> = \Theta(x^0-y^0)  \Big( D(x-y) - D(y-x) \Big). \tag{1}$$
In contrast, other sources (see e.g., this answer and references therein), define the retarded propagator as
$$G_R = \Theta(x^0-y^0) \left< \phi(x) \phi(y) \right> =  \Theta(x^0-y^0) D(x-y) \tag{2}$$
What motivates these two clearly different definitions and, in particular, the much more complicated definition in Eq. 1?

The propagator given in Eq. 2 makes perfect sense. It's the probability amplitude to find the particle at $x=(t_x,\vec x)$ if it starts at $y=(t_y,\vec y)$ and is only nonzero if $t_x>t_y$. So we only take into account how a particle propagates to a different location at a later moment in time.
The propagator in Eq. 1 is stranger. It also contains the amplitude described above. But then we add to this amplitude the amplitude that the particle was at $y=(t_y,\vec y)$ if is now at $x=(t_x,\vec x)$. (By using the Heaviside function we make sure that $\phi(x)$ generates a state at an earlier moment in time. Hence the second term in Eq. 1 $\propto D(y-x) $ is the amplitude that the particle was at $y$ if it is now at $x$.)
So the propagator in Eq. 2 is something we can immediately understand why the propagator in Eq. 1 is quite unintuitive. Why does it make sense to consider the propagator in Eq. 1 and what's the physical difference between the two?

REPLY [2 votes]: I'm not an expert in the subject, but my immediate thought regarding the object defined in Eq. 2 is that it is not Lorentz-invariant.
As shown by Tong and Peskin, the two-point correlation function $\langle \phi(y) \phi(x) \rangle$ does not actually vanish outside the light cone.  In essence, there can be correlations between the values of the field at two points $x$ and $y$ without those points being causally connected.
However, if $(x-y)^2<0$ then by appropriate choice of Lorentz transformation we can change the sign of $x^0-y^0$, meaning that $\theta(x^0-y^0)\langle \phi(y)\phi(x)\rangle$ would vanish in one frame and not in another.
Of course, the equations you cite are just definitions.  Eq. 2 is fine, it's just that anything physical will end up being written as some Lorentz-invariant combination of those objects.  With the definition given in Eq. 1, we need not worry about such things because $\langle[\phi(y),\phi(x)]\rangle = 0$ outside the light cone, meaning that the sign of $x^0-y^0$ cannot be changed by Lorentz transformation and so everyone will agree on whether the propagator vanishes or not.


TITLE: Can lasers cause glare or reflections?
QUESTION [1 upvotes]: When going to the movie theater I sometimes notice reflections from the projector on my glasses and I find it very annoying. What happens is that the light hits the back of my glasses and i can see a reflection of it.

I am gonna watch a movie which is going to be projected using imax lasers. Would laser projectors produce glare of a similar nature?
How do laser lights work? Do the photons only go straight? Then how is it visible as a beam from the side sometimes?

REPLY [0 votes]: Having a laser projector tells you very little about how the projector actually works. The only thing to know for sure is that the light source is some kind of laser (as opposed to e.g. Xenon arc lamps). The advantage is usually that they can be made a lot brighter and with a wider color gamut (they can generate more colors).
But, the light travels no straighter than light from a regular projector. Simply because under normal circumstances (no black holes in the vicinity, no gradients of refractive index, etc.), light always travels in a straight line.
What you see in your glasses is probably stray light from imperfections in the beam path and scattered light from the dust in the air. These things will unfortunately still be there with a laser projector.


TITLE: What are perfect squares in general rings?
QUESTION [3 upvotes]: The definition of a square element can be made for arbitrary rings $R$:

$p\in R$ is a square if there is an $a\in R$ with $p = aa$.

The definition of a perfect square element for arbitary rings – inspired by perfect square numbers – would presuppose the concept of an integral element in arbitrary rings:

$p\in R$ is a perfect square if there is an integral element $a\in R$ with $p = aa$.

Note, that for  $R = \mathbb{Z}$ being a square and being a perfect square are equivalent, while for $R = \mathbb{Q}$ it's not.
But the definition of integral elements:

$p\in R$ is integral over $A$, a subring of $R$, if there are $n \geq 1$ and $a_{j}\in A$ such that $p^{n}+a_{n-1}p^{n-1}+\cdots
 +a_{1}p+a_{0}=0$

depends on a subring $A \subset R$. So also the definition of perfect squares depends on a subring $A \subset R$.
My questions are:

Is there an "absolute" (not relative) definiton of being a perfect square (like there are absolute definitons of being a prime element or a square element)?
Does the concept of being a perfect square (with respect to a subring $A \subset R$) play an important role in ring theory/algebra?
Are there sensible reasons that the name "perfect square" had been chosen (instead of the more self-explanatory name "integral square")?

REPLY [5 votes]: A perfect square (in normal use) is an element of $\Bbb Z$ which is the square of an element of $\Bbb Z$. That is, it is a square in $\Bbb Z$. Even if you look in some larger ring, say $\Bbb Q$, then the perfect squares are precisely the elements of $\Bbb Z$ which are squares in $\Bbb Z$. The fact that $\Bbb Z$ is a subring of $\Bbb Q$ is crucial to the definition - indeed, the point of the term "perfect square" is that for $\Bbb Z \subset R$ - e.g. $R = \Bbb C, \Bbb R, \Bbb Q$ - then the set of "perfect squares" in $R$ is the same for any choice of $R$ - it's the set of squares in $\Bbb Z$.
So why we have the term "perfect square" at all? Fundamentally, it's because if you call them "square numbers" then you get people saying "but $2$ is a square number it's the square of $\sqrt 2$". "Integral square" would work, except that this a concept you teach to 8 year olds and "perfect", unlike "integral", is a word they already know.
But if you're at the point where you can say "a square in the ring $R$", then the term "perfect square" isn't adding any mathematical content: it's just shorthand for "a square in the ring $\Bbb Z$".


TITLE: Is there a hilbert curve equivalent for circles?
QUESTION [3 upvotes]: Is there a space-filling curve that has the same properties of a hilbert curve (two points close in 1D are close in 2D) but grows in a circular shape instead of a rectangular one?

REPLY [3 votes]: I am not aware of a such curve construction, but the next best thing is probably using a suitable mapping between square and circle disc and then transform the square based hilbert curve onto a circle.
Maybe something like this: Conformal mapping circle onto square (and back)


TITLE: Equation of heat conduction for spherical solid
QUESTION [0 upvotes]: What will be the correct Equation of heat conduction for a homogeneous spherical solid with constant thermal diffusivity K and no heat source?

REPLY [1 votes]: For a medium with constant thermal diffusivity $K$ and no heat source the heat equation is
$$
\nabla^2 T(\pmb r, t)=\frac{1}{K}\frac{\partial T(\pmb r, t)}{\partial t}
$$
and in spherical coordinates the equation becomes
$$
{1 \over r^2}{\partial \over \partial r}\!\left(r^2 {\partial T \over \partial r}\right)
\!+\!{1 \over r^2\!\sin\theta}{\partial \over \partial \theta}\!\left(\sin\theta {\partial T \over \partial \theta}\right)
\!+\!{1 \over r^2\!\sin^2\theta}{\partial^2 T \over \partial \phi^2}=\frac{1}{K}\frac{\partial T}{\partial t}
$$


TITLE: Prove that $f$ is holomorphic in $U$ (second part)
QUESTION [1 upvotes]: This is the natural continuation of this other question. I have this exercise

Let $L:=a+\Bbb R b$, for some $a,b\in\Bbb C$, a straight line in $\Bbb C$. Also let $f\in C(U,\Bbb C)$ is holomorphic on $U\setminus L$. Show that $f$ is holomorphic on all of $U$. (Hint: Morera's theorem.)

I want to know is the proof below is correct or if there is some mistake or it is incomplete.
The version of Morera's theorem that I used here use only triangular paths $\Delta$, where a triangular path is a closed path defined by the straight segments between three different points whenever the region bounded by the triangle belongs to the domain of the function.
It says

Morera's theorem: if $\int_\Delta f dz=0$ for every triangular path in $U$ then $f$ is analytic in $U$.

The proof:
Observe that we need to consider just the analyticity of $f$ on $U\cap L$.
Then we only need to see what happen to the paths that contain some point of $U\cap L$. There are three kinds of triangular paths $\Delta$ that contain a point of $U\cap L$:

$|\Delta\cap L|=1$, then necessarily this point is a vertex of $\Delta$.
$|\Delta\cap L|=2$, then there is a segment of $L$ inside the area bounded by $\Delta$.
$|\Delta\cap L|=\mathfrak c$, then there is a whole segment of $\Delta$ contained in $L$.

Let a homotopy $H\in C^1([0,1]^2,U)$ such that each $H([0,1],t)$ is triangular. Now I want to show that if $f$ is continuous then the map $t\mapsto\int_{\Delta_t}f(z)\, dz$ is also continuous, where $\Delta_t:=H([0,1],t)$. Then
$$
\left|\int_{\Delta_{t_1}} f(z)\, dz-\int_{\Delta_{t_2}}f(z)\, dz\right|=\left|\int_0^1(f\circ H)(s,t_1)\partial_1H(s,t_1)\, ds-\int_0^1(f\circ H)(s,t_2)\partial_1H(s,t_2)\, ds\right|\\
=\left|\int_0^1\big((f\circ H)(s,t_1)\partial_1H(s,t_1)-(f\circ H)(s,t_2)\partial_1H(s,t_2)\big)\, ds\right|\\
\le\max_{s\in[0,1]}|(f\circ H)(s,t_1)\partial_1H(s,t_1)-(f\circ H)(s,t_2)\partial_1H(s,t_2)|\tag1
$$
Now observe that the function $g:=(f\circ H)(\cdot,\cdot)\partial_1 H(\cdot,\cdot)$ is uniformly continuous (because it is a continuous function with compact domain) so for every chosen $\epsilon>0$ there is a $\delta>0$ such that
$$
|(s,t_1)-(s,t_2)|=|t_1-t_2|<\delta\implies |g(s,t_1)-g(s,t_2)|<\epsilon\tag2
$$
Consequently the map $t\mapsto\int_{\Delta_t}f(z)\, dz$ is continuous. Thus if we choose $H$ such that $\Delta_t\subset U\setminus L$ for all $t\in[0,1)$ and $\Delta_1\subset U$ then
$$
\int_{\Delta_1}f(z)\, dz=\lim_{t\to 1}\int_{\Delta_t}f(z)\, dz=\lim_{t\to 1}0=0\tag3
$$
Now observe that this hold for every triangular path on $U$ of kind $1$ and $3$ because the interior of the region bounded by these triangles belongs completely to $U\setminus L$, so we can define such $H$ in any case. By last observe that a triangle of kind $2$ can be decomposed in two or three triangular paths of kind $1$ and $3$.
Then $(3)$ holds for every $\Delta_1\subset U$ so, by Morera's theorem, $f$ is analytic in $U$.$\Box$

REPLY [1 votes]: Your proof is not correct. You only considered triangles such that one of its vertices belong to $L$. That's far from being the general case.