The Ordered Matrix Dirichlet for State-Space Models

Niklas Stoehr${}^{\text{{\tipaencoding Q}}}$ Benjamin J.
Radford${}^{\text{{\tipaencoding N}}}$ Ryan Cotterell${}^{\text{{\tipaencoding
Q}}}$ Aaron Schein${}^{\text{{\tipaencoding@}}}$ ${}^{\text{{\tipaencoding
Q}}}$ETH Zurich ${}^{\text{{\tipaencoding N}}}$UNC Charlotte
${}^{\text{{\tipaencoding@}}}$The University of Chicago
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Many dynamical systems in the real world are naturally described by latent
states with intrinsic ordering, such as “ally”, “neutral”, and “enemy”
relationships in international relations. These latent states manifest through
countries’ cooperative versus conflictual interactions over time. _State-space
models (SSMs)_ explicitly relate the dynamics of observed measurements to
transitions in latent states. For discrete data, SSMs commonly do so through a
state-to-action _emission matrix_ and a state-to-state _transition matrix_.
This paper introduces the Ordered Matrix Dirichlet (OMD) as a prior
distribution over ordered stochastic matrices wherein the discrete
distribution in the $k^{\textrm{th}}$ row is stochastically dominated by the
$(k\\!+\\!1)^{\textrm{th}}$, such that probability mass is shifted to the
right when moving down rows. We illustrate the OMD prior within two SSMs: a
hidden Markov model, and a novel dynamic Poisson Tucker decomposition model
tailored to international relations data. We find that models built on the OMD
recover interpretable ordered latent structure without forfeiting predictive
performance. We suggest future applications to other domains where models with
stochastic matrices are popular (e.g., topic modeling), and publish user-
friendly code.

Figure 1: Action types in international relations data are ordered along a
cooperation-to-conflict axis. Our model infers latent states with an ordering
that reflects the ordering in observed actions. _Emission matrix_ : more
conflictual actions (indexed by higher $a$), are generated by more conflictual
latent states (higher $k$). _Transition matrix_ : latent states transition to
neighboring ones. _CDFs_ : The Ordered Matrix Dirichlet enforces that the CDF
of the $k^{\textrm{th}}$ discrete distribution is always greater than the
$(k\\!+\\!1)^{\textrm{th}}$, so that probability mass in the stochastic matrix
shifts right when moving down rows.

##  1 INTRODUCTION

In many modeling settings and application domains, some aspect of the
observation space has intrinsic ordering. For example, in international
relations, observed interactions between countries can be ordered on a
conflict-to-cooperation axis, ranging from “provide aid” to “fight”
(Goldstein,, 1992; Schrodt,, 2008). This ordering should ideally be reflected
in any state space used to summarize or describe observed interactions. For
example, more conflictual actions like “fight” or “threaten” might be more
likely between countries in an “enemy” state than those in an “ally”state
(Schrodt,, 2006). In this example, the latent states represent relationship
statuses between countries, ordered from “ally” to “enemy”, and reflect the
conflict-to-cooperation ordering of the observed actions. We might expect
states to transition to other states over time in a way that also reflects
their intrinsic ordering. Allies rarely become enemies from one moment to the
next, but rather (de-)escalate gradually, passing first through intermediate
states (Davis and Stan, 1984a, ).

State-space models (SSMs) are statistical models that explicitly relate time-
varying measurements to latent states, such that patterns and trends in the
observed space are attributable to transitions between states over time. For
discrete data, the canonical form of such models is based on two stochastic
matrices: the _emission matrix_ , which describes how latent states generate
observations, and the _transition matrix_ , which describes how states
transition to other states over time. This general formulation does not
intrinsically promote any ordering of the latent states, whose indices are
arbitrary and subject to “label switching” (Richardson and Green,, 1997;
Stephens,, 2000).

To promote some sense of ordering in the state space, researchers sometimes
constrain the transition matrix to take only banded or “left-right-left” forms
(Schrodt,, 2006; Netzer et al.,, 2008; Randahl and Vegelius,, 2022), whereby
states only transition to adjacent states. This constraint substantially
restricts the expressiveness of the model while still not ensuring a well-
defined ordering of the latent states that reflects the ordering in the
observation space.

This paper introduces a novel prior distribution over stochastic matrices, the
Ordered Matrix Dirichlet (OMD), and demonstrates it as a key ingredient in
SSMs with well-ordered state spaces. An OMD random variable is a stochastic
matrix whose rows are discrete distributions that sum to 1, and whose
$k^{\textrm{th}}$ row is stochastically dominated by the
$(k\\!+\\!1)^{\textrm{th}}$, so that probability mass shifts to the right when
moving down the matrix (Fig. 1). We define the OMD distribution implicitly via
a stick-breaking construction that ensures the desired ordering property by
sorting Beta-distributed auxiliary variables. As we show, when the OMD is
selected as a prior over both emission and transition matrices in an SSM, the
inferred latent states have an intrinsic ordering that reflects ordering in
the observation space.

To demonstrate and evaluate the OMD as a prior, we construct and apply two
different SSMs—(1) a simple hidden Markov model (HMM) that we apply to
synthetic data where the ground-truth latent structure is known, and (2) a
novel dynamic version of Bayesian Poisson Tucker decomposition (Schein et al.,
2016b, ), which we apply to international relations data of country-to-country
interactions.

We compare each of these models to a baseline that differs only in its prior
over the transition and emission matrices—instead of the OMD, it makes the
standard assumption of rows being independently Dirichlet-distributed, which
we term the Standard Matrix Dirichlet (SMD). We find that the models based on
the OMD are more readily interpretable than those based on the SMD while still
performing comparably and sometimes better in forecasting and imputation
tasks. In synthetic experiments, the OMD model is much more effective at
recovering ground-truth latent structure, while on international relations
data, the OMD model exhibits superior forecasting performance over both the
SMD model and an additional baseline we introduce that constrains the
transition matrix to be banded.

After setting up notation and providing background on SSMs in § 2, we formally
introduce the OMD in § 3 and motivate it as a prior within SSMs. In § 4 we
discuss posterior inference for OMD models using Pyro (Bingham et al.,, 2018).
We then provide results from a suite of synthetic data experiments in § 5, and
present a case study on international relations data in § 6. Finally, we
discuss broader connections in § 7, and summarize our conclusions in § 8.

##  2 STATE-SPACE MODELS

State-space models (SSMs) describe the evolution of time-indexed measurements
$\boldsymbol{y}^{(t)}$ in terms of corresponding latent states
$\boldsymbol{\lambda}^{(t)}$ (Kalman,, 1960). SSMs assume that patterns and
trends in the observed measurements, typically only noisily realized, are
attributable to transitions between latent states.

##### Basic Form

This paper considers a subset of SSMs frequently used to model discrete or
non-negative data. Consider a non-negative vector-valued measurement
$\boldsymbol{y}^{(t)}\in\mathbb{R}_{+}^{A}$ at discrete time step $t$. Using
the international relations example in the introduction,
$\boldsymbol{y}^{(t)}$ might measure the counts of $A$ different action types
taken between some pair of countries during time step $t$. We use $a\in[A]$ to
index into this vector, so that $y^{(t)}_{a}$ is an entry, and refer to $a$ as
an _action_ or _action type_ throughout.

The SSMs we consider connect observed measurements to vector-valued, non-
negative latent states $\boldsymbol{\lambda}^{(t)}\in\mathbb{R}_{+}^{K}$ under
the following assumption:

$\displaystyle\mathbb{E}[y^{(t)}_{a}]$
$\displaystyle\propto\sum_{k=1}^{K}\lambda^{(t)}_{k}\underbrace{\phi_{ka}}_{\text{emission}}$
(1)

where $\phi_{ka}\in[0,1]$ is an entry in the discrete distribution
$\boldsymbol{\phi}_{k}$ which sums to one over actions
$\sum_{a=1}^{A}\phi_{ka}\\!=\\!1$ and is itself the $k^{\textrm{th}}$ row of
the state-to-action _emission matrix_ $\Phi$.

The states then evolve over time under the assumption

$\displaystyle{\mathbb{E}}[\lambda^{(t)}_{k}]$
$\displaystyle\propto\sum_{k^{\prime}=1}^{K}\lambda^{(t\\!-\\!1)}_{k^{\prime}}\underbrace{\pi_{k^{\prime}k}}_{\text{transition}}$
(2)

where $\pi_{k^{\prime}k}\in[0,1]$ is an entry in the discrete distribution
$\boldsymbol{\pi}_{k^{\prime}}$ which is itself the
$(k^{\prime})^{\textrm{th}}$ row of the state-to-state _transition matrix_
$\Pi\in[0,1]^{K\times K}$.

Many SSMs follow this basic form in 1 and 2, such as hidden Markov models
(HMMs), discrete dynamical systems (Schein et al., 2016a, ), or more complex
SSMs as presented in § 6.2. The key feature of these models is that they
involve an emission $\Phi$ and transition matrix $\Pi$ which are both
_(row-)stochastic matrices_.

##### What is a “State”?

There are differences in the SSM literature on whether the “state” at time
step $t$ _is_ the vector $\boldsymbol{\lambda}^{(t)}$, or whether each element
$\lambda^{(t)}_{k}$ of the vector describes the relevance of one of $k\in[K]$
“states”. These two interpretations coincide when $\boldsymbol{\lambda}^{(t)}$
is a one-hot vector, placing non-zero mass on only one element, as in HMMs.
However, these interpretations diverge in more general settings. We adopt both
senses of the word “state” in this paper, referring to
$\boldsymbol{\lambda}^{(t)}$ as the _complex “state” of the overall system_ at
$t$ but also understanding it as a _mixture over $K$ simple “states”_.

##### Dirichlet Priors

Researchers often place prior distributions over model parameters either as a
way to encode structural assumptions about the state space or to fit models
using Bayesian inference (or both). The conventional prior for row-stochastic
matrices assumes that rows are independently Dirichlet distributed, what we
will refer to as the Standard Matrix Dirichlet (SMD). A draw
$\boldsymbol{\phi}\sim\textrm{Dir}(\boldsymbol{\alpha})$ from a Dirichlet
distribution with concentration parameter
$\boldsymbol{\alpha}\in\mathbb{R}_{+}^{A}$ is a discrete distribution over $A$
categories, $\sum_{a=1}^{A}\phi_{a}\\!=\\!1$.

##### Banded Constraints

One commonly-used constraint is that the transition matrix is banded along its
diagonal, such that $\pi_{k^{\prime}k}=0$ if $|k-k^{\prime}|>b$ for some
bandwidth $b$ (often set to 1); see the middle plot of Fig. 2. This constraint
encodes the assumption that states only excite or transition to nearby states
at subsequent time steps. Such an assumption is motivated, for example, when
the desired state space represents ordered stages of escalation in
international conflict (Schrodt,, 2006; Randahl and Vegelius,, 2022). For
purposes of comparison, we introduce a prior distribution called the Banded
Matrix Dirichlet (BMD) that enforces this constraint (§ B.2).

Figure 2: Stochastic matrices sampled from the Standard Matrix Dirichlet
(SMD), Banded Matrix Dirichlet (BMD) and Ordered Matrix Dirichlet (OMD).
Neither the SMD nor the BMD adhere to the stochastic dominance property in Eq.
3, as evidenced by overlapping CDFs.

##  3 THE ORDERED MATRIX DIRICHLET

Many dynamical systems in the real world are naturally described by latent
states with intrinsic ordering, such as in international relations, where the
relationship status of two countries might escalate from “ally” to “enemy”
only gradually, first passing through intermediate states like “neutral”. In
addition to constraining how states transition over time, this ordering may
further reflect ordering in the observed actions between countries, with
countries in more conflictual latent states (e.g., “enemy”) taking more
conflictual actions towards each other (e.g., “fight”), and countries in more
cooperative states (e.g., “ally”) taking more cooperative actions (e.g.,
“provide aid”).

The state space in the basic model formulation given in Eq. 1 and Eq. 2 is not
intrinsically ordered. Specifically, the row indices $k\in[K]$ of the emission
and transition matrices are arbitrary and bear no intrinsic information. As
mentioned in the previous section, constraining the transition matrix to be
banded does impart some information to the index $k$. However, its
interpretation is circular: $k$ is some state that transitions to other states
$\\{k^{\prime}:|k-k^{\prime}|\leq b\\}$, whose interpretation is similarly
defined with respect to $k$. Moreover, while banding the transition matrix
promotes some ordering of latent states, it does not promote one that
necessarily reflects the ordering in observed actions.

To overcome these limitations, this section introduces a novel prior over the
transition and emission matrices that ensures an intrinsically well-ordered
state space, one that both reflects the ordering in observed actions and the
ordering in latent state transitions. We first operationalize our notion of
ordering in terms of stochastic dominance, and then construct a probability
distribution with support over the subset of stochastic matrices that obey
this notion.

###  3.1 Ordering by Stochastic Dominance

Considering first the emission matrix, each row $\boldsymbol{\phi}_{k}$
represents a discrete distribution over ordinal actions types. Intuitively, we
might say the rows are well-ordered if probability mass shifts to the right
when moving down rows, or equivalently, when the $k^{\textrm{th}}$
distribution places more weight on earlier action types than the
$(k\\!+\\!1)^{\textrm{th}}$. This intuition is formalized by the notion of
_stochastic dominance_ (Davidson,, 2017). Define the cumulative distribution
function (CDF) for the $k^{\textrm{th}}$ discrete distribution to be
$\mathrm{CDF}_{k}(a)\triangleq\sum_{a^{\prime}=1}^{a}\phi_{ka^{\prime}}$.
Then, the $k^{\textrm{th}}$ distribution is stochastically dominated by the
$(k\\!+\\!1)^{\textrm{th}}$ if

$\mathrm{CDF}_{k}(a)\geq\mathrm{CDF}_{k\\!+\\!1}(a)\,\,\textrm{ for all }a$
(3)

We refer to a stochastic matrix as _well-ordered_ if Eq. 3 holds for all rows
$k$. For the $K\times A$ emission matrix, this means higher rows place more
mass on earlier actions types. For the $K\times K$ transition matrix, this
notion encapsulates many but not all banded structures (see Fig. 2), while
further allowing for much more flexible “down-and-right” transition shapes
(see Fig. 4).

Figure 3: _Heatmaps:_ Averages of 10 samples from the Ordered Matrix Dirichlet
with varying concentration ${\boldsymbol{\alpha}}\in\mathbb{R}_{+}^{A}$.
“Column mean” refers to
$\overline{\phi_{a}}=\frac{1}{K}\sum_{k=1}^{K}{\boldsymbol{\phi}}_{ka}$ and
shows that probability mass is asymmetric to the right. _Simplex plots:_ Each
point is a discrete distribution over $A=3$ classes. The $K=3$ points ( , , )
connected by a line represent one sample from the OMD. We observe ordered
transitions from the lower left $[1,0,0]$ to the top $[0,0,1]$ corner of the
simplex. In Fig. 12 in App. C, we present unordered sample trajectories from
the SMD in comparison.

###  3.2 Ordered Matrix Dirichlet (OMD) Distribution

We now introduce a new probability distribution that has support over only the
matrices described above. The OMD distribution is defined by two parameters,
_concentration_ ${\boldsymbol{\alpha}}\in\mathbb{R}_{+}^{A}$ and _height_ $K$.
An OMD random variable $\Phi\sim\textrm{\text{OMD}}(K,{\boldsymbol{\alpha}})$
is a $K\times A$ matrix that is row-stochastic and well-ordered, as shown in
Proposition 3.1.

We define the OMD distribution implicitly via Algorithm 1, which generates OMD
variates. This algorithm builds on the stick-breaking construction of the
standard Dirichlet distribution (Gelman et al.,, 2013, p. 583). A Dirichlet
random variable $\phi\sim\textrm{Dir}(\boldsymbol{\alpha})$ can be generated
iteratively, one entry at a time, via Beta auxiliary variables. First, draw
$\phi_{1}\sim\textrm{Beta}(\alpha_{1},\sum_{a>1}\alpha_{a})$. Then for
$a=2,\dots,A\\!-\\!1$ draw
$\beta_{a}\sim\textrm{Beta}(\alpha_{a},\sum_{a^{\prime}>a}\alpha_{a^{\prime}})$
and set
$\phi_{a}\leftarrow\beta_{a}\,(1-\sum_{a^{\prime}<a}\phi_{a^{\prime}})$.
Finally, set $\phi_{A}\leftarrow 1-\sum_{a^{\prime}<A}\phi_{a^{\prime}}$.
Intuitively, this construction iteratively “breaks” off some amount of
remaining probability mass (the “stick”), where the Beta variables determine
the size of the breaks.

Algorithm 1 iteratively constructs $K$ discrete distributions over $A$
categories using the same basic idea. For each category $a$ in succession, it
samples $K$ Beta variables (lines 3 and 8) to determine the size of the
“breaks” in the $K$ “sticks”. However, it further sorts the Beta variables
(lines 5 and 10), so that the largest “break” of the remaining “stick” is
always taken by the first “stick” ($k=1$), the second largest is always taken
by the second “stick” ($k=2$), and so on. In so doing, it generates a well-
ordered stochastic matrix, as stated below.

######  Proposition 3.1 (OMD random variables are well-ordered).

The OMD has support over only row-stochastic matrices that obey the ordering
property given beneath Eq. 3, such that for any two rows $k<k^{\prime}$ and
any $a$

$\sum_{a^{\prime}\leq a}\phi_{ka^{\prime}}\geq\sum_{a^{\prime}\leq
a}\phi_{k^{\prime}a^{\prime}}$ (4)

Proof: See § B.1.

##### Lack of Analytic Form

We define the OMD implicitly by construction and do not (yet) know any
analytic form for its probability density function (PDF), which involves
integrating over products of Beta order statistics. We leave further
investigation into the OMD’s PDF, moments, and other analytic properties for
the future. As we show in the next section though, its lack of analytic form
does not hamper posterior inference with modern probabilistic programming.

##### What is “Dirichlet” about the OMD?

The OMD’s name reflects its definition as a minimal modification to the stick-
breaking construction of the Standard Matrix Dirichlet—if we remove the blue
lines (5 and 10), then Algorithm 1 corresponds exactly to the SMD, which
simply generates $K$ independent Dirichlet variates (with no ordering). The
name reflects this alone—it is not the case (to our knowledge) that the $K$
discrete distributions, which are dependent under the OMD via the sort
operation, are marginally or conditionally Dirichlet distributed.

Algorithm 1 Ordered Matrix Dirichlet

1:Input: height $K$, concentration
${\boldsymbol{\alpha}}\in\mathbb{R}_{+}^{A}$

2:for $k=1,\,\dots\,,K$ do

3:
$\tilde{{\boldsymbol{\phi}}}_{k1}\sim\textrm{Beta}\left(\alpha_{1},\sum_{a=2}^{A}\alpha_{a}\right)$

4:end for

5:$(\phi_{11},\dots,\phi_{K1})\leftarrow\textsc{Sort}\Big{(}(\tilde{\phi}_{11},\dots,\tilde{\phi}_{K1})\Big{)}$

6:for $a=2,\,\dots\,,\,A-1$ do

7: for $k=1,\,\dots\,,K$ do

8:
$\tilde{\beta}_{ka}\sim\textrm{Beta}\left(\alpha_{a},\sum_{a^{\prime}=a+1}^{A}\alpha_{a^{\prime}}\right)$

9: end for

10:
$(\beta_{1a},\dots,\beta_{Ka})\leftarrow\textsc{Sort}\Big{(}(\tilde{\beta}_{1a},\dots,\tilde{\beta}_{Ka})\Big{)}$

11: for $k=1,\,\dots\,,K$ do

12:
${\boldsymbol{\phi}}_{ka}\leftarrow\left(1-\sum_{a^{\prime}=1}^{a-1}\phi_{ka^{\prime}}\right)\beta_{ka}$

13: end for

14:end for

15:for $k=1,\,\dots\,,K$ do

16: ${\boldsymbol{\phi}}_{kA}\leftarrow
1-\sum_{a^{\prime}=1}^{A-1}\phi_{ka^{\prime}}$

17:end for

18:Output: OMD variate $\Phi\in\mathbb{R}_{+}^{K\times A}$

##### Concentration Parameter

The OMD is parameterized by its concentration
${\boldsymbol{\alpha}}\in\mathbb{R}_{+}^{A}$. Fig. 3 visualizes OMD samples
$\Phi$ for different settings of ${\boldsymbol{\alpha}}$. For symmetric
${\boldsymbol{\alpha}}=(\alpha_{0},\dots,\alpha_{0})$, one might expect
samples $\Phi$ to distribute probability mass across the matrix evenly.
However, we observe otherwise, that mass often skews to the right,
particularly for larger $\alpha_{0}$; we speculate this relates to the sorting
operation. Although unappealing, this does not mean the OMD is inherently
asymmetric, but rather that non-trivial settings of ${\boldsymbol{\alpha}}$
may be required to promote symmetry in the prior. Despite this, in practice,
we find that samples from the posterior are often symmetric.

##### Label Switching

The problem of “label switching” (Stephens,, 2000; Murphy,, 2012, p. 841)
arises in (ad)mixture models when the indices $k$ of latent states are
arbitrary, such that permuting them gives the same joint probability under the
model. This issue can hamper interpretation and prevents one from averaging
parameters across posterior samples without first aligning the states (e.g.,
using the Hungarian matching algorithm (Kuhn,, 1955)). The models we have
discussed, which place an OMD prior over the emission matrix, have
intrinsically ordered states that are not prone to label switching (see Fig.
11 in App. C). Although we do not view this as the main benefit of the OMD, it
is a welcome side effect that facilitates easier interpretation and permits
direct averaging of posterior parameters without post-hoc, potentially error-
prone alignment methods.

##  4 MCMC INFERENCE WITH PYRO

The OMD integrates nicely with modern probabilistic programming frameworks
like Pyro (Bingham et al.,, 2018; Phan et al.,, 2019).111We open-source our
code with tutorials and examples at https://github.com/niklasstoehr/ordered-
matrix-dirichlet Although we do not have an analytic form for its PDF, we are
able to build and perform efficient gradient-based MCMC on a range of OMD-
based models by implementing Algorithm 1. We regard the OMD as a modeling
motif that blends white- and black-box approaches in a way that was only
recently made feasible by advances in scientific computing.

We use Pyro’s implementation of the No-U-Turn Sampler (NUTS, Hoffman and
Gelman,, 2014), a variant of Hamiltonian Monte Carlo (HMC, Duane et al.,,
1987), to perform approximate posterior inference in OMD-based models. As with
any MCMC method, this returns a set of $S$ posterior samples of model
parameters $\\{\Pi^{(s)},\Phi^{(s)},\dots\\}_{s=1}^{S}$ which collectively
approximate the posterior distribution. In practice, we take $S=$$1000$
samples after $200$ burn-in samples.

NUTS relies on first-order gradient information of the model’s unnormalized
log joint density. Our implementation in Pyro takes gradients of the OMD
density implicitly via backpropagation through the stick-breaking
construction. Although the sort operation is not fully differentiable, it is
piece-wise linear and sub-differentiable (Boyd and Vandenberghe,, 2004;
Blondel et al.,, 2020; Tim Vieira,, 2021). We can view sort as a combination
of two operations: first the non-differentiable argsort obtains permutation
indices, then the differentiable gather applies the permutation. In the
backward pass, the permutation of indices is simply reversed to match their
original positions, obviating the need to differentiate through argsort. For
this reason, the sorting of Beta variates in the construction of the OMD does
not hinder gradient-based MCMC methods for inference.

Figure 4: SMD versus OMD at forecasting. _Column 1_ : Stylized ground-truth
transition matrices. _Columns 2-3_ : OMD recovers the transition matrices
while the SMD suffers from label switching. _Columns 4:_ OMD is better at
forecasting than SMD but worse at imputation. _Columns 5:_ OMD recovers the
latent states while the SMD suffers from label switching.

##  5 SYNTHETIC DATA EXPERIMENTS

We conduct experiments with synthetic data to better understand and evaluate
the behavior of SSMs with OMD priors. In particular, we generate datasets
using hidden Markov models (HMMs) with roughly diagonal emission matrices and
a range of stylized transition matrices—i.e., “banded”, “bonbon” and
“triangle”, all displayed in the left column of Fig. 4. The “bonbon”, for
instance, represents a realistic scenario for political event data, where
“neutral” states fluctuate but “ally” and “enemy” states are nearly absorbing.
To each of these datasets, we fit an HMM with OMD priors and compare its
performance to a baseline HMM with SMD priors.

We generate multiple datasets for each parameter setting using $10$ random
seeds, where each dataset comprises $N=$ $10000$ sequences of length $T=10$,
and where a single observation takes one of $A=10$ ordinal values. We further
consider two settings: one with all $N=$ $10000$ sequences and a few-shot
setting where models are fit to only $N=100$. We also generate random train-
test splits to evaluate two different forms of prediction:

  1. (1)

Imputation: We mask a random $30$% of all observations which models impute
during inference.

  2. (2)

Forecasting: We designate the first $70$% of time steps for training and the
latter $30$% for testing. Models are fit to the training set, then used to
forecast the test observations.

##### Qualitative Results.

We first compare how well the two models recover known ground-truth latent
structure. As expected, the OMD model reliably recovers the shape of the true
transition matrix while the SMD model does not, often exhibiting label
switching; see the first 3 columns of Fig. 4 for examples. As a simple
quantitative measure of this, we can also calculate the mean absolute error
(MAE) between the true latent states and the inferred ones. The last column of
Fig. 4 reports the error on forecasting future latent states where the OMD
model is substantially better; this is unsurprising and simply confirms that
the OMD’s states are well-ordered while the SMD’s are label-switched.

##### Predictive Results.

The $4^{\textrm{th}}$ column of Fig. 4 reports MAE on forecasting future
observations. The OMD model performs at least as well as the SMD model in all
settings, and sometimes substantially better, as when the true transition
matrix is banded ($1^{\textrm{st}}$ row). By contrast, the imputation results
in Fig. 14 in App. C show the OMD model performing substantially worse than
the SMD model in most settings. We speculate that the strong inductive bias
imparted by the OMD prior helpfully regularizes the model’s forecasts while
overly restricting its imputation ability. Intriguingly, the one setting where
the OMD model has superior imputation performance is when the true transition
matrix is banded, which accords with the forecasting results.

Figure 5: ICEWS event data showing interactions between Armenia and Azerbaijan
over monthly time steps. Figure 6: Imputation and forecasting evaluation on
held-out ICEWS data, over $10$ runs with random seed. We fit the DPT model
with different parametrizations (SMD, BMD, OMD) of the emission $\Phi$ and
transition $\Pi$ matrix. We find that OMD does not significantly reduce
predictive results suggesting that the imposed constraints fit the given data.
Figure 7: Posterior mean of parameters of Dynamic Poisson Tucker model fitted
to ICEWS subset of 2020. (A) The latent matrix $\Psi$ indicates country-
community activity. Armenia and Azerbaijan are standing out as they are
involved mostly in one community (B) Interactions between latent communities.
We find that communities $c=1$, $c=3$ and $c=4$ predominantly interact between
themselves. (C) Selected community interactions over time: $c_{1}=4\rightarrow
c_{2}=4$ are in conflict, while $c_{1}=2\rightarrow c_{2}=5$ are mostly
neutral. (D) Latent emission matrix $\Phi$ representing global state-to-action
probabilities. Thanks to the ordering, we know that state $k=3$ represents
conflictual relationships. (E) Latent transition matrix $\Pi$ representing
smooth state-to-state transition probabilities.

##  6 CASE STUDY: POLITICAL EVENTS

In this section, we give a case study on building an OMD-based SSM to analyze
international relations event data.

###  6.1 ICEWS Political Events Data

We consider political event data from the Integrated Crisis Early Warning
System (ICEWS) dataset (Boschee et al.,, 2015). ICEWS event data comprise
millions of micro-records of the form “country $i$ took action $a$ to country
$j$ at time $t$” that are machine-extracted from digital news archives. The
country actors $i$ and $j$ and action types $a$ are coded to follow the
Conflict and Mediation Event Observations (CAMEO) ontology (Schrodt,, 2012).

##### Ordered Actions

CAMEO specifies $20$ high-level action types, depicted in Fig. 5, that are
naturally ordered on a conflictual-to-cooperative axis—specifically, they are
each assigned a value on the expert-elicited Goldstein scale (Goldstein,,
1992), where the most cooperative action, “provide aid”, has a value of
$+7.0$, and the most conflictual action, “use unconventional mass violence”,
has a value of $-10.0$.

##### 4-mode Count Tensor

Following Schein et al., 2016b , we represent the data as a count tensor
${\boldsymbol{Y}}\in\mathbb{N}_{0}^{V\times V\times A\times T}$, where an
element $y^{(t)}_{i\xrightarrow{a}j}$ is the number of times country $i$ took
action $a$ to country $j$ during time step $t$. We consider $V=249$ countries,
$A=20$ action types (ordered by Goldstein values), and $T=72$ months. The
$A\times T$ slice of this tensor corresponding to all interactions between
$i\equiv\textsc{Armenia}$ and $j\equiv\textsc{Azerbaijan}$ is visualized in
Fig. 5.

###  6.2 Dynamic Poisson Tucker model

Our model assumes each count $y^{(t)}_{i\xrightarrow{a}j}$ is Poisson
distributed:

$\displaystyle y^{(t)}_{i\xrightarrow{a}j}$
$\displaystyle\sim\textrm{Pois}\left(\delta_{a}\,\delta^{(t)}\sum_{k=1}^{K}\lambda^{(t)}_{i\xrightarrow{k}j}\underbrace{\phi_{ka}}_{\text{emission}}\right)$
(5)

where $\phi_{ka}$ is an entry in the state-to-action emission matrix, the
parameters $\delta_{a}$ and $\delta^{(t)}$ are action- and time-scaling
coefficients, and $\lambda^{(t)}_{i\xrightarrow{k}j}$ represents how well the
$k^{\textrm{th}}$ state describes the relationship $(i\rightarrow j)$ at time
$t$. Eq. 5 conforms to the basic form given in Eq. 1, where here the
measurements and states are specifically tensor-valued.

Our model further assumes that $\lambda^{(t)}_{i\xrightarrow{k}j}$ decomposes
so that

$\displaystyle\sum_{k=1}^{K}\lambda^{(t)}_{i\xrightarrow{k}j}\phi_{ka}\equiv\sum_{c_{1}=1}^{C}\psi_{c_{1}i}\sum_{c_{2}=1}^{C}\psi_{c_{2}j}\sum_{k=1}^{K}\lambda^{(t)}_{c_{1}\xrightarrow{k}c_{2}}\phi_{ka}$
(6)

where $\psi_{c_{1}i}$ and $\psi_{c_{2}j}$ represent the rate at which
countries $i$ and $j$ participate in _latent communities_ $c_{1}$ and $c_{2}$,
respectively, and $\lambda^{(t)}_{c_{1}\xrightarrow{k}c_{2}}$ then represents
how well the $k^{\textrm{th}}$ state describes the inter-community
relationship $(c_{1}\rightarrow c_{2})$ at time $t$. The multilinear form in
Eq. 6 corresponds to a Tucker decomposition Tucker, (1964), where the
$\lambda^{(t)}_{c_{1}\xrightarrow{k}c_{2}}$ values collectively form the core
tensor $\boldsymbol{\Lambda}^{(t)}\in\mathbb{R}_{+}^{C\times C\times K}$ at
time $t$. In this setting, the core tensor can also be interpreted as a
tensor-valued _state (of the whole system)_.

We then model the evolution of the core tensor over time as

$\displaystyle\lambda^{(t)}_{c_{1}\xrightarrow{k}c_{2}}$
$\displaystyle\sim\textrm{Gam}\left(\tau_{0}\,\sum_{k^{\prime}=1}^{K}\lambda^{(t-1)}_{c_{1}\xrightarrow{k^{\prime}}c_{2}}\underbrace{\pi_{k^{\prime}k}}_{\text{transition}},\,\tau_{0}\right)$
(7)

which follows the form of Poisson–Gamma Dynamical Systems (Schein et al.,
2016a, ) while conforming to Eq. 2.

We place non-informative gamma priors over the parameters
$\delta_{a},\,\psi_{c_{1}i},\,\psi_{c_{2}j}\stackrel{{\scriptstyle\textrm{iid}}}{{\sim}}\textrm{Gam}(\alpha_{0},\alpha_{0})$,
a dynamic prior over
$\delta^{(t)}\sim\textrm{Gam}(\tau_{0}\delta^{(t\\!-\\!1)},\,\tau_{0})$, and
set $\tau_{0}=\alpha_{0}=1$.

Finally, with all of the aforementioned structure the same, we then consider
three different settings for the priors over the transition $\Pi$ and emission
$\Phi$ matrices: (1) OMD+OMD, where both are drawn from the OMD, (2) SMD+SMD,
where both are drawn from the SMD, and (3) SMD+BMD, where $\Phi$ is drawn from
the SMD and $\Pi$ is drawn from the Banded Matrix Dirichlet (BMD), as defined
in § B.2.

###  6.3 Experiments and Results

To further understand and evaluate the OMD we fit the three above-mentioned
versions of the Dynamic Poisson Tucker (DPT) model to ICEWS data and compare
their qualitative and predictive performance. We use the same hyperparameters
for all models with $C=5$ and $K=3$.

##### Predictive Evaluation

Following the design in § 5, we create $10$ train-test splits that randomly
mask observations for imputation and withhold later time steps for
forecasting. In addition to MAE, we evaluate performance using _scaled
pointwise predictive density (SPPD)_ , a measure between $0$ and $1$ where
higher is better, which we define in § B.3.

Fig. 6 reports the imputation and forecasting results for each of the three
models. As in the synthetic experiments, we see that the OMD model is better
than the SMD at forecasting but worse at imputation. Similarly, the BMD model
is also better than the SMD at forecasting but worse at imputation. This
strengthens our belief that the OMD’s inductive bias regularizes its forecasts
while overly restricting its imputation ability, since the BMD exhibits the
same pattern, and their two inductive biases are similar. That being said, the
OMD is much more flexible than the BMD, which may explain why it outperforms
the BMD in both forecasting and imputation.

##### Qualitative Exploration

To qualitatively inspect its inferred latent structure, we fit the OMD model
to the fully-observed dataset. Fig. 7 visualizes the _posterior mean_ of
inferred model parameters for the time period of 2020. Since the model is not
prone to label switching, we can inspect the posterior mean, as opposed to
inspecting single (often arbitrary) sample. Fig. 7A visualizes the country-
community matrix $\Psi$. We observe that Armenia and Azerbaijan are
predominantly involved in community $c=4$. By visualizing a slice of the core
tensor in Fig. 7B, we see that this community mostly interacts with itself. By
visualizing in Fig. 7C the slice $\lambda^{(t)}_{4\xrightarrow{k}4}$, we see
which states $k$ best describe community $c=4$’s self-interactions over time.
In mid 2020, the most active state is $k=3$. We immediately know it represents
a conflictual relationship since its index $k$ is high. This is confirmed by
the emission matrix in Fig. 7D where we see that state $k=3$ places most of
its mass on “fight”. Finally, we visualize the transition matrix in Fig. 7E
and find that, unfortunately, transitioning out of state $k=3$ seems unlikely.

We also visualize inferred latent structure from a DPT model with $K=6$ and
$C=20$ fitted to ICEWS data from a longer time range 2015–2020 and present the
results in Fig. 9.

##  7 DISCUSSION

##### International Relations

As alluded to throughout, this work was largely motivated by datasets,
modeling approaches, and core concepts in the field of international relations
(IR). The notion of _escalation_ —that countries only gradually transition to
conflict through an orderly sequence of intermediate states—is fundamental to
how scholars organize and understand political events (Davis and Stan, 1984b,
). Theoretical accounts for why countries fight attempt to characterize a
sequence of intermediate states that rational actors would transition through
on their way to war (Snyder,, 1984; Fearon,, 1995; Jervis,, 2017). A similar
perspective underlies empirical approaches. The earliest attempts to digitize
international affairs into “event data” were explicitly couched in the
framework of escalation—the very first sentence of Azar, (1980) reads: “As
students of politics and political science, we should and we do care about the
events which lead to war…”

The principal challenge in the data-intensive study of international relations
is the inherent sparsity and missingness of event data, which provide only a
scattered glimpse at the underlying structures we seek to reason about. This
paper follows an empirical tradition of encoding strong inductive biases into
statistical models of event data which encourage their inferred structure to
accord with theoretical notions, like “escalation” (Schrodt,, 2006; Anders,,
2020; Randahl and Vegelius,, 2022). While much of the previous work focuses on
constraining (specifically, banding) the transition structure between “states”
to encourage orderly dynamics, the key idea in this paper is to draw further
on the ordinal nature of observed action types. There is a steadily-growing
literature on models for dyadic event data that has made exciting advances
while still mostly treating action types as unordered (O’Connor et al.,, 2013;
Schein et al.,, 2015; Minhas et al.,, 2016). In parallel, there has been
recent work on inferring latent intensity scales Terechshenko, (2020); Stoehr
et al., (2022) that imbue actions with a richer or more data-driven sense of
ordering. We are eager for these threads to continue to cross, as they have in
this work.

##### Other Models and Other Domains

The Ordered Matrix Dirichlet as a modeling motif is applicable beyond
international relations and SSMs. We include in App. C a brief exploration of
other OMD-based models we have built, with illustrative results on other
datasets. Building these models in Pyro is easy, often only requiring a few
lines of code, which facilitates this exploration. Fig. 8 summarizes four
different models, all of which (and more) are available in the code we have
open-sourced; we describe them here too.

  1. (1)

We build an ordered form of Poisson–Gamma Dynamical Systems (PGDS) (Schein et
al., 2016a, ) placing an OMD prior over the transition and emission matrices.
PGDS was originally introduced to model ICEWS data, but treats actions as
unordered.

  2. (2)

We use an HMM with OMD-distributed emission and transition matrices to model
the observed global change of temperature. In this model, noisy temperature
changes are related to ordered latent states indicative of “warming” and
“cooling” periods that transition gradually.

  3. (3)

Even simpler, we experiment with a Markov chain model consisting of a single
state-to-state transition matrix to model sleep cycles. Sleep cycles typically
transition step-by-step from wake (W) to rapid eye movement (REM) stages (Pan
et al.,, 2012).

  4. (4)

We modify Latent Dirichlet Allocation (LDA) (Blei et al.,, 2003) to place an
OMD prior over the topic-word matrix. While word types are canonically viewed
as unordered, we imbue them with ordering by sorting them on “semantic axes”
(An et al.,, 2018), for instance from negative to positive words. The model
then infers ordered topics that reflect this semantic axis, similar to the
model of Stoehr et al., (2023).

We can imagine many more applications that motivate well-ordered state-space
models, like modeling product life cycles (Arvidsson,, 2019) or customer-
company relationships (Netzer et al.,, 2008). Beyond SSMs, admixture models,
like LDA, are fundamentally based on stochastic matrices and used in
population genetics (Pritchard et al.,, 2000), stochastic block models
(Airoldi et al.,, 2008), recommender systems (Gopalan et al.,, 2015), among
many other areas.

##  8 CONCLUSION

This paper introduced the Ordered Matrix Dirichlet (OMD) distribution as a
prior distribution over well-ordered stochastic matrices in state-space models
(SSMs). Models built on the OMD have intrinsically ordered states that reflect
ordering in the observed data. These models are more readily interpretable and
usable, as they are not prone to label switching, while still being
competitive on predictive tasks. The OMD integrates nicely with modern
probabilistic programming frameworks, making it easy to build and fit OMD-
based models. While this paper’s motivation is rooted in the concepts and data
of international relations, the motifs presented here have broad applicability
to domains with ordinal data and models based on stochastic matrices.

#### Acknowledgments

We would like to thank Kevin Du and the anonymous reviewers for valuable
feedback on the manuscript and Tim Vieira for his input on order statistics
and sorting. Niklas Stoehr is supported by the Swiss Data Science Center
(SDSC).

## References

  * Airoldi et al., (2008) Airoldi, E. M., Blei, D., Fienberg, S., and Xing, E. (2008).  Mixed membership stochastic blockmodels.  Advances in neural information processing systems, 21.
  * An et al., (2018) An, J., Kwak, H., and Ahn, Y.-Y. (2018).  SemAxis: A lightweight framework to characterize domain-specific word semantics beyond sentiment.  In Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics, volume 1.
  * Anders, (2020) Anders, T. (2020).  Territorial control in civil wars: Theory and measurement using machine learning.  Journal of Peace Research, 57(6):701–714.
  * Arvidsson, (2019) Arvidsson, R. (2019).  On the use of ordinal scoring scales in social life cycle assessment.  The International Journal of Life Cycle Assessment, 24(3):604–606.
  * Azar, (1980) Azar, E. E. (1980).  The conflict and peace data bank (copdab) project.  Journal of Conflict Resolution, 24(1):143–152.
  * Bingham et al., (2018) Bingham, E., Chen, J. P., Jankowiak, M., Obermeyer, F., Pradhan, N., Karaletsos, T., Singh, R., Szerlip, P., Horsfall, P., and Goodman, N. D. (2018).  Pyro: Deep universal probabilistic programming.  Journal of Machine Learning Research.
  * Blei et al., (2003) Blei, D. M., Ng, A. Y., and Jordan, M. I. (2003).  Latent Dirichlet allocation.  Journal of Machine Learning Research, 3:993–1022.
  * Blondel et al., (2020) Blondel, M., Teboul, O., Berthet, Q., and Djolonga, J. (2020).  Fast differentiable sorting and ranking.  International Conference on Machine Learning.
  * Boschee et al., (2015) Boschee, E., Lautenschlager, J., O’Brien, S., Shellman, S., Starz, J., and Ward, M. (2015).  ICEWS Coded Event Data.
  * Boyd and Vandenberghe, (2004) Boyd, S. and Vandenberghe, L. (2004).  Convex optimization.  Cambridge University Press.
  * Davidson, (2017) Davidson, R. (2017).  Stochastic dominance.  In The New Palgrave Dictionary of Economics, pages 1–7.
  * (12) Davis, P. K. and Stan, P. (1984a).  Concepts and Models of Escalation.  RAND Corporation.
  * (13) Davis, P. K. and Stan, P. J. E. (1984b).  Concepts and models of escalation: A report from the Rand Strategy Assessment Center.  Rand.
  * Duane et al., (1987) Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987).  Hybrid Monte Carlo.  Physics Letters B, 195(2):216–222.
  * Fearon, (1995) Fearon, J. D. (1995).  Rationalist explanations for war.  International organization, 49(3):379–414.
  * Gelman et al., (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian data analysis.  Chapman and Hall/CRC.
  * Gelman et al., (2014) Gelman, A., Hwang, J., and Vehtari, A. (2014).  Understanding predictive information criteria for Bayesian models.  Statistics and Computing, 24(6):997–1016.
  * Goldstein, (1992) Goldstein, J. (1992).  A conflict-cooperation scale for WEIS events data.  The Journal of Conflict Resolution, 36(2):369–385.
  * Gopalan et al., (2015) Gopalan, P., Hofman, J. M., and Blei, D. M. (2015).  Scalable recommendation with hierarchical Poisson factorization.  In UAI, pages 326–335.
  * Hoffman and Gelman, (2014) Hoffman, M. D. and Gelman, A. (2014).  The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15(1):1593–1623.
  * Jervis, (2017) Jervis, R. (2017).  Perception and misperception in international politics: New edition.  Princeton University Press.
  * Kalman, (1960) Kalman, R. E. (1960).  A new approach to linear filtering and prediction problems.  Journal of Basic Engineering, 82(1):35–45.
  * Kuhn, (1955) Kuhn, H. W. (1955).  The Hungarian method for the assignment problem.  Naval Research Logistics Quarterly, 2(1-2):83–97.
  * Minhas et al., (2016) Minhas, S., Hoff, P. D., and Ward, M. D. (2016).  A new approach to analyzing coevolving longitudinal networks in international relations.  Journal of Peace Research, 53(3):491–505.
  * Murphy, (2012) Murphy, K. P. (2012).  Machine learning: A probabilistic perspective.  MIT Press, Cambridge, MA.
  * Netzer et al., (2008) Netzer, O., Lattin, J. M., and Srinivasan, V. (2008).  A hidden Markov model of customer relationship dynamics.  Marketing Science, 27(2):185–204.
  * O’Connor et al., (2013) O’Connor, B., Stewart, B., and Smith, N. (2013).  Learning to extract international relations from political context.  In Proceedings of the 51st Annual Meeting of the Association for Computational Linguistics, volume 1, pages 1094–1104.
  * Pan et al., (2012) Pan, S.-T., Kuo, C.-E., Zeng, J.-H., and Liang, S.-F. (2012).  A transition-constrained discrete hidden Markov model for automatic sleep staging.  BioMedical Engineering OnLine, 11(1):52.
  * Phan et al., (2019) Phan, D., Pradhan, N., and Jankowiak, M. (2019).  Composable effects for flexible and accelerated probabilistic programming in NumPyro.  arXiv, 1912.11554.
  * Pritchard et al., (2000) Pritchard, J. K., Stephens, M., and Donnelly, P. (2000).  Inference of population structure using multilocus genotype data.  Genetics, 155(2):945–959.
  * Randahl and Vegelius, (2022) Randahl, D. and Vegelius, J. (2022).  Predicting escalating and de-escalating violence in Africa using Markov models.  International Interactions, pages 1–17.
  * Richardson and Green, (1997) Richardson, S. and Green, P. J. (1997).  On Bayesian analysis of mixtures with an unknown number of components.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59(4):731–792.
  * Schein et al., (2015) Schein, A., Paisley, J., Blei, D. M., and Wallach, H. (2015).  Bayesian Poisson tensor factorization for inferring multilateral relations from sparse dyadic event counts.  In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1045–1054.
  * (34) Schein, A., Wallach, H., and Zhou, M. (2016a).  Poisson-Gamma dynamical systems.  In Advances in Neural Information Processing Systems, volume 29.
  * (35) Schein, A., Zhou, M., Blei, D. M., and Wallach, H. (2016b).  Bayesian Poisson Tucker decomposition for learning the structure of international relations.  In Proceedings of the 33rd International Conference on Machine Learning, volume 48, pages 2810–2819.
  * Schrodt, (2008) Schrodt, P. (2008).  Kansas event data system (KEDS).
  * Schrodt, (2012) Schrodt, P. (2012).  CAMEO: Conflict and mediation event observations event and actor codebook.  Parus Analytics.
  * Schrodt, (2006) Schrodt, P. A. (2006).  Forecasting conflict in the Balkans using hidden Markov models.  In Programming for Peace, pages 161–184. Springer Netherlands.
  * Snyder, (1984) Snyder, G. H. (1984).  The security dilemma in alliance politics.  World politics, 36(4):461–495.
  * Stephens, (2000) Stephens, M. (2000).  Dealing with label switching in mixture models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 62(4):795–809.
  * Stoehr et al., (2023) Stoehr, N., Cotterell, R., and Schein, A. (2023).  Sentiment as an ordinal latent variable.  In European Chapter of the ACL (EACL).
  * Stoehr et al., (2022) Stoehr, N., Hennigen, L. T., Valvoda, J., West, R., Cotterell, R., and Schein, A. (2022).  An ordinal latent variable model of conflict intensity.  In arXiv, volume 2210.03971.
  * Terechshenko, (2020) Terechshenko, Z. (2020).  Hot under the collar: A latent measure of interstate hostility.  Journal of Peace Research, 57(6):764–776.
  * Tim Vieira, (2021) Tim Vieira (2021).  On the distribution function of order statistics.
  * Tucker, (1964) Tucker, L. R. (1964).  The extension of factor analysis to three-dimensional matrices.  In Gulliksen, H. and Frederiksen, N., editors, Contributions to mathematical psychology., pages 110–127. Holt, Rinehart and Winston.

##  Appendix A IMPACT STATEMENT

We emphasize that our models are intended for research purposes and empirical
insights. They should not be blindly deployed for automated decision-making
processes. The used ICEWS data may contain biases that are potentially
reinforced by our modeling assumptions. The experiments with real-world event
data in § 6.3 were conducted on an NVIDIA TITAN RTX GPU. The experiments with
synthetically generated data in § 5 can be run on a local M1 CPU with $64$ GB
of RAM in less than $10$ minutes. Limiting factors are the selected
hyperparameter sizes for the latent states $K$ and communities $C$, as well as
the number of time series $N$ and their length $T$. We discuss further model
limitations in § 8 and § 3.

##  Appendix B SUPPLEMENTARY TECHNICAL DETAILS

###  B.1 Proof of Proposition 3.1

######  Proposition B.1 (OMD random variables are well-ordered).

The OMD has support over only row-stochastic matrices that obey the ordering
property given beneath Eq. 3, such that for any two rows $k<k^{\prime}$ and
any $a$

$\sum_{a^{\prime}\leq a}\phi_{ka^{\prime}}\geq\sum_{a^{\prime}\leq
a}\phi_{k^{\prime}a^{\prime}}$ (8)

Proof: For $a=1$, $\phi_{k1}>\phi_{k^{\prime}1}$ is true by construction (line
5 of Algorithm 1). For $a=2$, by the definition in line 12,
$\phi_{k2}=(1-\phi_{k1})\,\beta_{k2}$, and therefore the CDF at $a=2$ equals
$\phi_{k1}-\phi_{k1}\beta_{k2}+\beta_{k2}$. It suffices to show that
$\phi_{k1}-\phi_{k1}\beta_{k2}+\beta_{k2}>\phi_{k^{\prime}1}-\phi_{k^{\prime}1}\beta_{k^{\prime}2}+\beta_{k^{\prime}2}$,
since the remaining $a>2$ then follow by induction. Re-arranging terms,
$(\phi_{k1}\\!-\\!\phi_{k^{\prime}1})+(\beta_{k2}\\!-\\!\beta_{k^{\prime}2})>(\phi_{k1}\beta_{k2}-\phi_{k^{\prime}1}\beta_{k^{\prime}2})$,
which follows since we know by construction that
$\beta_{k2}>\beta_{k^{\prime}2}$ (line 10), and all terms
$\phi_{k1},\phi_{k^{\prime}1},\beta_{k2},\beta_{k^{\prime}2}$ are between 0
and 1.

###  B.2 Details of the Banded Matrix Dirichlet (BMD)

In this section, we elaborate on the Banded Matrix Dirichlet (BMD). For
simplicity, we consider a square matrix $\Pi\in[0,1]^{K\times K}$, but the BMD
can be non-square as well. We assume that the $k^{\textrm{th}}$ state can only
be excited by its directly neighboring states, $(k-1)^{\textrm{th}}$ and
$(k+1)^{\textrm{th}}$, as well as by itself (Schrodt,, 2006; Randahl and
Vegelius,, 2022). This results in a matrix whose non-zero elements are banded
along the diagonal following:

$\displaystyle\pi_{kk^{\prime}}=\begin{cases}\pi_{k}^{{(\nearrow)}}&\textrm{
if }k^{\prime}=k+1\,\,\textrm{ (escalating)}\\\
\pi_{k}^{{(\searrow)}}&\textrm{ if }k^{\prime}=k-1\,\,\textrm{
(descalating)}\\\ \pi_{k}^{{(\circ)}}&\textrm{ if }k^{\prime}=k\,\,\textrm{
(steady)}\\\ 0&\textrm{ otherwise }\end{cases}$ (9)

Finally, we place a Dirichlet prior over the three non-zero elements in each
$k^{\textrm{th}}$ row

$\displaystyle(\pi_{k}^{{(\nearrow)}},\pi_{k}^{{(\searrow)}},\pi_{k}^{{(\circ)}})$
$\displaystyle\sim\textrm{Dir}(\alpha_{0}^{{(\nearrow)}},\alpha_{0}^{{(\searrow)}},\alpha_{0}^{{(\circ)}})$
(10)

Moreover, we can consider a wider bandwidth $b\geq 1$ so that components
$k^{\prime}\in\\{k-b,\dots,k+b\\}$ all excite $k$. An example of the full
vector might then look like

$\displaystyle{\boldsymbol{\pi}}_{k}=(0,\dots,0,\pi_{k}^{{(\searrow)}},\pi_{k}^{{(\circ)}},\pi_{k}^{{(\nearrow)}},0,\dots,0)$
(11)

###  B.3 Details on the Scaled Pointwise Predictive Density

_Scaled pointwise predictive density (SPPD)_ is defined as

$\textrm{SPPD}=\exp\Big{(}\tfrac{1}{|\mathcal{I}|}\sum_{\mathbf{i}\in\mathcal{I}}\log\Big{[}\tfrac{1}{S}\sum_{s=1}^{S}\textrm{Pois}\big{(}y_{\mathbf{i}};\,\mu_{\mathbf{i}}^{(s)}\big{)}\Big{]}\Big{)}$
(12)

where $\mathbf{i}$ is the multi-index of an entry $y_{\mathbf{i}}$ in the
tensor—e.g., $\mathbf{i}=(i,j,a,t)$—and $\mathcal{I}$ is the set multi-indices
corresponding to all entries in the test set. The term
$\mu_{\mathbf{i}}^{(s)}$ is the Poisson rate in Eq. 5 as given by the
$s^{\textrm{th}}$ posterior sample of model parameters. SPPD is the same as
LPPD (Gelman et al.,, 2014), but scaled by $\tfrac{1}{|\mathcal{I}|}$ and
exponentiated so it is always between 0 and 1, where higher is better.

###  B.4 Relevant Links

Code accompanying this paper
https://github.com/niklasstoehr/ordered-matrix-dirichlet
Integrated Crisis Early Warning System (ICEWS)
https://dataverse.harvard.edu/dataverse/icews
Goldstein Scale
https://parusanalytics.com/eventdata/cameo.dir/CAMEO.SCALE.txt
Conflict and Mediation Event Observations (CAMEO)
https://parusanalytics.com/eventdata/cameo.dir/CAMEO.09b6.pdf

##  Appendix C SUPPLEMENTARY PLOTS

Figure 8: Different models with Dirichlet-sampled latent matrices fitted on
data exhibiting ordinal dynamics. The Latent Dirichlet Allocation (LDA) has no
temporal dimension, but similarly comprises a stochastic matrix describing
word distributions per latent topic. If we order the observed vocabulary of
words by the words’ sentiment score, the Ordered Matrix Dirichlet (OMD) can
recover topics representative of sentiment levels. In all settings, we find
that the OMD yields more easily interpretable stochastic matrices than the
Standard Matrix Dirichlet (SMD). Figure 9: Posterior mean of parameters of
Dynamic Poisson Tucker model, with $K=6$ latent states and $C=20$ latent
communities, fitted to full temporal range (2015-2020) of ICEWS data. We find
that the probability mass of the transition matrix is centered along the
diagonal revealing step-wise (de-)escalatory dynamics. There is high
probability of staying in state $k=6$ indicating that conflictual
relationships may be hard to escape. The country-community affiliation matrix
$\Psi$ provides no information on whether communities represent allies or
enemies per se. To obtain this information, we interact the country-community
matrix with the core tensor
$\psi^{(\rightarrow)}_{c_{1}i}\sum_{c2=1}^{C}\sum_{j=1}^{V}\psi^{(\leftarrow)}_{c_{2}j}\lambda^{(t)}_{c_{1}\xrightarrow{k}c_{2}}$
for specific choice of $k$ and $t$.

Figure 10: Descriptive statistics showing total number of interactions between
countries in ICEWS data from 2015 to 2020. The rows and columns are sorted by
the total number of actions a country is involved in. Note that we omit self-
targeted actions as indicated by the blank diagonal.

Figure 11: Recovering ground truth structures in transition and emission
matrices of a state-space model. Conventionally, rows are samples
independently from a (standard) Dirichlet distribution. This can result in
label switching making the latent states (topics) difficult to interpret. This
is particularly problematic if states are ordinal, e.g., representing “ally”,
“neutral” and “enemy” relations.

Figure 12: Samples from the Standard Matrix Dirichlet (SMD). Each point in the
triangle plot represents a sample from a Dirichlet over $A=3$ classes. The
$K=3$ points connected by a line represent an (unordered) sample from the SMD.

Figure 13: Imputation results of synthetic data experiments. As discussed in §
5, we generate time series with different ground truth transition structures:
“banded”, “bonbon”, “triangle”. We fit a Hidden Markov Model (HMM) to a train
set of these data and evaluate imputation performance on a test set. In
contrast to the forecasting experiments (Fig. 4), SMD \+ SMD outperforms OMD
\+ OMD in two out of three cases on observed states. In contrast to
forecasting, imputation does not necessarily require a model with temporal
dynamics and the ordered transition matrix does not help. As expected, OMD \+
OMD performs better at imputing latent states because it circumvents label
switching.

| action type
---
a
| action
---
name
| Goldstein
---
value
0 | provide aid | 7.0
1 | engage material cooperation | 6.0
2 | yield | 5.0
3 | express intent cooperate | 4.0
4 | engage diplomatic cooperation | 3.5
5 | appeal | 3.0
6 | consult | 1.0
7 | make public statement | 0.0
9 | investigate | -2.0
10 | disapprove | -2.0
11 | reject | -4.0
12 | reduce relations | -4.0
13 | demand | -5.0
14 | threaten | -6.0
15 | protest | -6.5
16 | coerce | -7.0
17 | exhibit force posture | -7.2
18 | assault | -9.0
19 | fight | -10.0
20 | unconventional mass violence | -10.0
Figure 14: Ordered CAMEO action types with assigned Goldstein values. We order
action types by Goldstein value first and, in case of a tie, by CAMEO ID
second.


# Modeling the 3-micron Class Er-Doped Fluoride Fiber Laser with a Cubic
Energy Transfer Rate Dependence

William Bisson, Alexandre Michaud, Pascal Paradis, Réal Vallée, and Martin
Bernier

###### Abstract

We propose an energy transfer model with a cubic atomic population dependence
to accurately model the behavior of various reported high-power erbium-doped
fluoride fiber lasers operating near 2.8 microns. We first show that the
previously introduced weakly interacting (WI) and strongly interacting (SI)
models are not adequate for precisely modeling such high-power erbium-doped
fluoride fiber lasers. We compare results obtained with the WI and SI models
to the proposed model by simulating 4 different highly doped (7 mol.%) fiber
lasers previously reported in the literature. Laser efficiencies and powers
are reproduced with great accuracy. In addition, four other independent fiber
laser systems based on erbium doping concentrations varying from 1-6 mol.% are
also simulated with good accuracy using the proposed model with the exact same
set of spectroscopic parameters, which confirms its validity for various
erbium doping concentrations. Redshifting of laser wavelength is also taken
into account by considering the full cross section spectra and computing
signal powers over several wavelength channels.

###### Index Terms:

Erbium-doped fiber laser, erbium, Er3+, laser theory, energy transfer, mid-
infrared, 2.8 $\mu$m, numerical modeling.

This work has been submitted to the IEEE for possible publication. Copyright
may be transferred without notice, after which this version may no longer be
accessible.

##  I Introduction

High-power erbium-doped fluoride fiber lasers operating near 2.8 microns are
undergoing rapid development due to their various high-end applications, for
instance in biomedicine, material processing and remote sensing [1]. A
commonly used pumping band for these lasers is the erbium transition from the
fundamental level ${}^{4}I_{15/2}$ to level ${}^{4}I_{11/2}$, which is around
976 nm [3, 5, 4, 6, 7, 8, 2]. The Stokes efficiency of erbium lasers pumped in
this band is around 33%, but it was found that colasing at 1.6 $\mu m$ can
help increasing this efficiency to 50% [2]. The Stokes efficiency was also
exceeded by using highly doped systems, which led to higher energy transfer
upconversion rates, and an optical-to-optical efficiency of 35.4% was reached
by Faucher et al. [3]. More recently, the record CW power of 42 W, with an
overall efficiency of 22.9%, was reported by Aydin et al. [4] by pumping from
both ends of the laser cavity and by improving the fiber laser assembly.

The modeling of erbium-doped fluoride fiber lasers operating around 2.8
microns has been a subject of interest for multiple years [10, 9, 12, 11], but
no satisfactory model was found to fit the multiple reported lasers with the
same set of spectroscopic parameters. Having a precise model is important
since it can guide experiments and namely help scaling up output powers. It
can also give access to underlying information, like atomic populations for
example [13]. A certain flaw actually exists in the literature regarding the
energy transfer rates, for which two sets of parameters were proposed. These
sets are named the weakly interacting (WI) and the strongly interacting (SI)
models [9]. The SI model relies on transfer rates that were derived from bulk
samples whereas the WI model uses the same general equations but with scaled
down values for the transfer rates. In addition, these sets of transfer rates
change with the total erbium concentration [10], which means they have to be
measured every time a new concentration is used. The WI model was first
introduced by Gorjan et al. [11] to obtain better agreement with experimental
values, while Li and Jackson used the same set of energy transfer parameters
to simulate multiple lasers [9]. Although the WI model provides a better
agreement with the experiment than does the SI model, it is however unable to
reproduce the experimental lasing efficiency properly at all pumping powers.
This can be seen in figure 1 where we show simulation results of the 20 W
laser from ref. [3] using the WI and SI models. The simulations in ref. [9]
were also unable to accurately reproduce the efficiency of this laser at all
pump powers.

Figure 1: Simulation of the 7 mol.% 2825 nm laser from ref. [3] with the WI
and SI models and a pump wavelength of 976 nm.

Both the SI and WI models rely on the same quadratic dependence of the ETU
processes as originally proposed by Grant [14]. It assumes infinitely fast
energy migration [15] and has been shown to be inappropriate in many
experimental cases [16].

In this paper, we propose an alternative handling of the energy transfer
mechanism for the accurate modeling of high-power diode pumped erbium-doped
fluoride fiber lasers operating around 2.8 microns. This approach is based on
previously reported models [17, 18, 15, 19] and combines the extreme cases of
energy transfer without the assistance of migration and energy transfer that
is purely assisted by migration. We show that the proposed energy transfer
model, which has a cubic dependence on excited atomic level population,
provides better agreement with experimental data than the WI and SI models for
4 different heavily-doped (7 mol.%) fiber lasers reported in the literature
[3, 4, 5, 6]. We show that such model can also be used, with the exact same
set of spectroscopic parameters, for modeling a 6 mol.% laser [7], an
unpublished 5 mol.% laser and two colasing (at 1.6 and 2.8 microns) systems
based on 1.5 and 1 mol.% erbium-doped fluoride fibers respectively [8, 2]. The
proposed model reproduces the efficiency drop in accordance to pump power
which is generally observed in heavily erbium-doped fluoride fiber lasers with
good accuracy, and is thus suitable for predictions related to power scaling
of such laser systems. It does not consider the clustering of ions, which
still has not been directly proven for erbium-doped fluoride fiber lasers
[20]. We believe the new insight provided by this alternative modeling could
pave the way to the accurate design of fiber laser systems operating over the
100W-level at 2.8 microns.

##  II Theoretical Model

###  II-A 1.6 microns and 2.8 microns erbium model

The energy level diagram containing the relevant radiative transitions and
energy transfer processes for the proposed model is shown in figure 2.

Figure 2: Energy level diagram and lifetimes of erbium-doped fluoride glass.
GSA = Ground State Absorption, ESA = Excited State Absorption, ETU = Energy
Transfer Upconversion, CR = Cross Relaxation, M = Migration.

The atomic populations are calculated using the following set of rate
equations [9, 12]:

$\displaystyle\derivative{N_{6}}{t}$
$\displaystyle=-N_{6}\tau_{6}^{-1}+R_{esa1}+f_{etu2},$ (1)
$\displaystyle\derivative{N_{5}}{t}$
$\displaystyle=\beta_{65}N_{6}\tau_{6}^{-1}-N_{5}\tau_{5}^{-1}-f_{cr},$ (2)
$\displaystyle\derivative{N_{4}}{t}$
$\displaystyle=\sum_{j=5}^{6}\beta_{j,4}N_{j}\tau_{j}^{-1}-N_{4}\tau_{4}^{-1},$
(3) $\displaystyle\derivative{N_{3}}{t}$
$\displaystyle=\sum_{j=4}^{6}\beta_{j,3}N_{j}\tau_{j}^{-1}-N_{3}\tau_{3}^{-1}+R_{esa2}$
$\displaystyle+f_{cr}+f_{etu1},$ (4) $\displaystyle\derivative{N_{2}}{t}$
$\displaystyle=\sum_{j=3}^{6}\beta_{j,2}N_{j}\tau_{j}^{-1}-N_{2}\tau_{2}^{-1}+R_{gsa}-R_{esa1}$
$\displaystyle-R_{se28}-2f_{etu2},$ (5) $\displaystyle\derivative{N_{1}}{t}$
$\displaystyle=\sum_{j=2}^{6}\beta_{j,1}N_{j}\tau_{j}^{-1}-N_{1}\tau_{1}^{-1}-R_{se16}+R_{se28}$
$\displaystyle-R_{esa2}+f_{cr}-2f_{etu1},$ (6) $\displaystyle N_{tot}$
$\displaystyle=\sum_{j=0}^{6}N_{i},$ (7)

which is solved in the steady-state regime by setting
$\derivative{N_{i}}{t}=0$. The atomic populations are then time independent
$(N_{i}(z))$. $\tau_{i}$ are the lifetimes and $\beta_{ij}$ are the branching
ratios. Their values for level 1 to level 5 were taken from [9], and values
for $\tau_{6}$ and $\beta_{6j}$ were taken from [12].

The pump transition rates around 976 nm are given by:

$\displaystyle R_{gsa}$
$\displaystyle=\sum_{\lambda}\xi\quantity[N_{0}\sigma_{02}-N_{2}\sigma_{20}][P_{+}+P_{-}]\Delta\lambda,$
(8) $\displaystyle R_{esa1}$
$\displaystyle=\sum_{\lambda}\xi\quantity[N_{2}\sigma_{26}-N_{6}\sigma_{62}][P_{+}+P_{-}]\Delta\lambda.$
(9)

The laser signal transition rates around 2.8 $\mu$m and 1.6 $\mu$m are given
by:

$\displaystyle R_{se16}$
$\displaystyle=\sum_{\lambda}\xi\quantity[N_{1}\sigma_{10}-N_{0}\sigma_{01}][P_{+}+P_{-}]\Delta\lambda,$
(10) $\displaystyle R_{esa2}$
$\displaystyle=\sum_{\lambda}\xi\quantity[N_{1}\sigma_{13}-N_{3}\sigma_{31}][P_{+}+P_{-}]\Delta\lambda,$
(11) $\displaystyle R_{se28}$
$\displaystyle=\sum_{\lambda}\xi\quantity[N_{2}\sigma_{21}-N_{1}\sigma_{12}][P_{+}+P_{-}]\Delta\lambda,$
(12)

Where $se$ stands for stimulated emission. The sum is over all spectral
channels, $\Delta\lambda$ is the spectral length of each channel and $P_{+}$
and $P_{-}$ are the power spectral densities in W/m propagating in the forward
and backward directions in the fiber respectively. $P_{+}$ and $P_{-}$ are
separated in an array containing a power density value for each wavelength
channel. $\sigma_{ij}$ are the wavelength dependent absorption and emission
cross sections. The $\xi$ factor is also wavelength dependent and is given by:

$\xi=\frac{\lambda\Gamma}{hcA_{\text{core}}},$ (13)

where $\lambda$ is the wavelength of the considered channel, $A_{\text{core}}$
is the fiber core area and $\Gamma$ is the power-filling factor. For pump
signals, which propagate in the cladding, the power-filling factor is given by
the square of the ratio of the core area to the cladding area of the fiber:

$\Gamma=\frac{A_{core}^{2}}{A_{clad}^{2}}=\frac{r_{core}^{2}}{r_{clad}^{2}}.$
(14)

For signal propagating in the core of the fiber, the power filling factor is:

$\Gamma=1-\frac{u^{2}}{V^{2}}\quantity(1-\frac{K_{\nu}^{2}(w)}{K_{\nu+1}(w)K_{\nu-1}(w)}),$
(15)

where $u$, $w$ are the core and cladding transverse parameters, respectively,
$V$ is the normalized frequency of the fiber and $K_{\nu}$ is the modified
Bessel function of the second kind of order $\nu$. In our simulations of
colasing systems, two wavelength ranges, spanning from 1.4 $\mu$m to 1.7
$\mu$m and from 2.6 $\mu$m to 3.0 $\mu$m, were used. For lasers that only emit
around $2.8~{}\mu$m, only the range from 2.6 $\mu$m to 3.0 $\mu$m was used,
which means the optical transitions around 1.6 $\mu$m were not considered for
these lasers. These power channels were separated in increments of
$\Delta\lambda=0.5$ nm. Using many channels provides the ability to let the
system converge by itself to a lasing wavelength, instead of setting it as an
input parameter [13].

The pump absorption cross sections $\sigma_{02}$ and $\sigma_{26}$ were taken
from [21] and their respective emission cross sections ($\sigma_{20}$ and
$\sigma_{62}$) were calculated using McCumber theory [22, 23] by using the
energy values from [24]. For the laser signal cross sections, $\sigma_{10}$
and $\sigma_{01}$ were taken from [25], $\sigma_{13}$ and $\sigma_{31}$ from
[12], while $\sigma_{21}$ and $\sigma_{12}$ were taken from [26]. The full
cross section spectra were used for the simulations.

According to the currently accepted energy transfer model, the three transfer
rates are given by:

$\displaystyle f_{etu1}^{G}$ $\displaystyle=W_{11}N_{1}^{2},$ (16)
$\displaystyle f_{etu2}^{G}$ $\displaystyle=W_{22}N_{2}^{2},$ (17)
$\displaystyle f_{cr}^{G}$ $\displaystyle=W_{05}N_{0}N_{5},$ (18)

Where $W_{11}$, $W_{22}$ and $W_{05}$ are constants given in ref. [9]. The
$N^{2}$ dependence was derived by Grant [14] considering infinitely fast
energy migration [15].

The power propagation equation is given by:

$\derivative{P_{\pm}}{z}=(g-\alpha)P_{\pm}$ (19)

where $P_{\pm}$ represent the power spectral densities in W/m propagating in
forward or backward directions, $g$ and $\alpha$ are the gain and background
loss at a given wavelength channel respectively and $dz$ is the length of each
fiber section, which was set to 4 cm or less depending on the simulated
system. The gain is given by:

$g=\Gamma\sum_{\text{rad}}(\sigma_{ji}N_{j}-\sigma_{ij}N_{i}),$ (20)

where the sum is on all radiative transitions. Spontaneous emission was
considered by adding a single photon in each wavelength channel [27].

###  II-B Energy Transfer Model

According to Dexter [17], the energy transfer rate between two ions is given
by the following equation [15]:

$W_{\text{DD/DA}}=\frac{C_{ik,j\ell}}{d^{6}}.$ (21)

Equation (21) assumes electric dipole-dipole interactions, which is the
dominant process in the case of mid-infrared photonics [19, 18]. $d$ is the
dipole distance and $C_{ik,j\ell}$ the Förster microparameter given by [15]:

$C_{ik,j\ell}=\frac{3c}{8\pi^{4}n^{2}}I_{ik,j\ell}.$ (22)

$I_{ik,j\ell}$ is the overlap integral of the process and is given by:

$I_{ik,j\ell}=\int\sigma^{e}_{ij}(\lambda)\sigma^{a}_{k\ell}(\lambda)d\lambda,$
(23)

where $\sigma^{e}_{ij}$ is the cross section of the emitting ion going from
level $i$ to $j$ and $\sigma^{a}_{k\ell}$ is the cross section of the
absorbing ion going from level $k$ to $\ell$. Equation (21) to (23) are valid
for any energy transfer process, and overlap integrals for ETU1, ETU2, CR and
migration processes are shown in table LABEL:table:overlap. The value for the
CR process was guessed by comparing energy transfer values with the WI model
since the $\sigma_{53}$ cross section is unknown.

TABLE I: Overlap integrals associated with each transfer process used in the theoretical model and references for the cross sections used.  Process | $I_{ik,j\ell}$ | Overlap [10-58m5] | Ref
---|---|---|---
ETU1 | $I_{10,13}$ | 1.0 | [12, 25]
ETU2 | $I_{20,26}$ | 3.1 | [21]
M1 | $I_{10,01}$ | 101 | [25]
M2 | $I_{20,02}$ | 4.5 | [21]
CR | $I_{53,01}$ | 0.7 | Guessed
M5 | $I_{50,05}$ | 2.8 | [26]

According to Zubenko [16], Grant’s energy transfer model [14] does not apply
in many experimental cases, since it assumes an infinitely fast energy
migration process. As can be seen in figure 1, we were unable to reproduce the
efficiency drop which was observed in ref. [3] with a quadratic dependence of
the energy transfer rate. We believe this points out the fact that Grant’s
model does not apply to the 3-microns class erbium-doped fluoride fiber laser
system. Other excited level dependencies were thus tested by simulating
multiple high-power fluoride fiber lasers found in the literature [3, 4, 5, 6,
7, 8, 2]. Simulation results for all of these lasers are shown in the next
section. It was found empirically that the following energy transfer rates for
ETU1, ETU2 and CR gave the most accurate simulation results:

$\displaystyle f_{etu1}$ $\displaystyle=a^{-6}C_{etu1}N_{1}^{3}$
$\displaystyle+\sqrt{\frac{\pi}{6}}a^{-6}\sqrt{C_{etu1}C_{m1}}N_{0}^{1/4}N_{1}^{11/4},$
(24) $\displaystyle f_{etu2}$ $\displaystyle=a^{-6}C_{etu2}N_{2}^{3}$
$\displaystyle+\sqrt{\frac{\pi}{6}}a^{-6}\sqrt{C_{etu2}C_{m2}}N_{0}^{1/4}N_{2}^{11/4},$
(25) $\displaystyle f_{cr}$ $\displaystyle=a^{-6}C_{cr}(N_{0}N_{5})^{3/2}$
$\displaystyle+\sqrt{\frac{\pi}{6}}a^{-6}\sqrt{C_{cr}C_{m5}}(N_{0}N_{5})^{3/2},$
(26)

where $C_{etu}$, $C_{cr}$ and $C_{m}$ are the Förster microparameters for ETU,
CR and migration respectively and $a=\quantity(\frac{4\pi}{3})^{-1/3}$. It was
found that adding a correction factor of 3% to $a$ gave the best results.
Equation (24) to (26) can be applied in the concentration range of at least 1
mol.% to 7 mol.%, which will be shown in the next section. The cubic
dependence of the energy transfer rate was found by assuming an uniform
distribution of ions [28]. The first term was inferred directly from Dexter’s
theory [17] and considers the donor to acceptor transfer process without the
assistance of migration [19]. The second term was inferred by using the
hopping regime of purely migration assisted energy transfer and was calculated
from the transfer rate given in ref. [18]. Combining these two terms to obtain
the total transfer rates was found to be necessary in order to obtain accurate
simulation results.

Note that for the hopping regime of migration assisted energy transfer to
apply for a specific energy transfer process, the condition
$\frac{C_{etu,cr}}{C_{m}}<<1$ must be verified [16]. Using equation (22), we
get:

$\frac{C_{etu,cr}}{C_{m}}=\frac{I_{ik,j\ell}^{\text{etu,cr}}}{I_{i^{\prime}k^{\prime},j^{\prime}\ell^{\prime}}^{\text{m}}}$
(27)

Using the overlap integral values in table LABEL:table:overlap, we get a ratio
of 0.01 for ETU1, 0.7 for ETU2 and 0.25 for CR, which justifies the use of the
hopping regime to compute the migration assisted transfer rate. For the second
term of the ETU transfer rates, the $N_{0}$ dependence comes from the donor
concentration dependence of energy migration transfer rate [15, 16] and was
found to be necessary in order to obtain accurate simulation results.

In the next section, we show that the proposed energy transfer model provides
improved modeling results as compared to the currently accepted WI and SI
models, especially at high pump powers in heavily-doped fibers, by simulating
8 different laser systems, all with the exact same set of spectroscopic
parameters.

##  III Results and Discussion

###  III-A 7 mol.% erbium-doped fiber lasers

The results obtained for the proposed model with 4 different 7 mol.% erbium-
doped fiber lasers are shown on figures 4 to 7. Unless otherwise specified,
for the spectrum of the background loss, we used data from LeVerreFluoré [29]
for a singlemode fluoride fiber, which is shown in figure 3. For each modelled
laser, this spectrum was normalized to fit the measured loss value at 2.8
$\mu$m, since they use different batches of fluoride fibers. Full spectra of
the fiber Bragg gratings (FBGs) were used if available. Otherwise, Gaussian
spectra were assumed.

Figure 3: Single mode fluoride fiber attenuation spectrum taken from
LeVerreFluoré’s website. A value of 35 dB/km was assumed for the pump
wavelengths. Figure 4: Simulation of the laser from [3] using the proposed
model and the WI and SI parameters. The intracavity splice loss and
attenuation at signal wavelength were fixed to 6.5% and 0.1 dB/m respectively.
The pump wavelength was fixed to 976 nm and the laser wavelength was 2825 nm.

Figure 4 shows the experimental signal power at 2825 nm vs. the absorbed pump
power at 976 nm for the laser in ref. [3] and simulations with the three
tested sets of energy transfer parameters. With the same set of cavity
parameters, the proposed model is much more accurate than the WI and SI
models. In fact, we were able to reproduce the efficiency change from 35% to
28% with good accuracy, which was not possible with the WI and SI parameters.
This efficiency change has not been reproduced that precisely before [9]. An
intracavity splice loss of 6.5%, which is between the average value of 3.5%
and the maximum value of 9.2% [3], gave the best results.

Figure 5: Simulation of the laser from [5] using the proposed model and the WI
and SI parameters. The attenuation at the signal wavelength was fixed to 0.17
dB/m. The pump wavelength was fixed to 976 nm, and the laser wavelength was
around 2824 nm because we applied a shift to the peak Bragg wavelength as
stated in ref. [5]. Loss induced by the FBG was also considered.

The model was then tested on the 5 W laser from ref. [5]. The results are
shown in figure 5. This laser has been properly simulated before, but only by
taking thermal effects into account [9], which was not necessary in our case.

Figure 6: Simulation of the laser from [4] using the proposed model and the WI
and SI parameters. The attenuation at the signal wavelength was fixed to 0.1
dB/m. The pump wavelength was changed from 969 nm to 981 nm according to the
measured values [4], and the laser wavelength was 2825 nm.

Results for the 42 W bidirectionally pumped fiber laser from ref. [4] are
shown in figure 6. For this laser, the pump wavelength shifted from 969 nm at
low power to 981 nm at the highest power. This wavelength shift was considered
in the simulation, which means different pump GSA and ESA values are used for
each pump power. The simulation with the proposed model gives an efficiency
that is too high up to about 125 W of launched pump power. We believe this
could be due to differences between the values used for the pump GSA and ESA
and their actual values for the fiber used in [4], which can be caused by
variation in the glass composition for example [21]. It could also simply be
caused by experimental error in these measured cross section spectra, but
their error are unknown [21]. Another cause could be error in the experimental
laser data or fiber properties in [4]. Nonetheless, the predicted laser powers
with the proposed model are much more accurate than the WI and SI models,
especially at high pump powers as shown in figure 6.

Figure 7: Simulation of the laser from [6] using the proposed model and the WI
and SI parameters. The attenuation at the laser wavelength was fixed to 0.2
dB/m. The pump wavelength was changed linearly from 970 nm at 1 W to 981 nm at
200 W according to the measured values, and the signal wavelength was 2938 nm.
Fusion splice loss for splices S2, S3 and S4 were fixed to 7%, 3% and 7%
respectively.

Results for the fourth simulated high-power 7 mol.% laser are shown in figure
7. The shift of the pump wavelength with power was also considered for this
laser and it is believed to be the main cause of the positive efficiency shift
with pump power of this laser [6]. It was necessary to modify the cross
section spectra around 2938 nm to reproduce the experimental results, which is
justified since this value is in the tail of the cross section which means it
has an higher uncertainty than values around 2825 nm for example [26]. The
$\sigma_{21}$ value taken from [26] was changed from $4.2\times 10^{-26}$ m2
to $3.0\times 10^{-26}$ m2 at 2938 nm, and $\sigma_{12}$ from $1.2\times
10^{-26}$ m2 to $1.5\times 10^{-26}$ m2 at the same wavelength.

###  III-B Other erbium concentrations

Another feature of the proposed model is that it can be used at other erbium
concentrations based on the same set of energy transfer rate equations. This
is evidenced in this section, where we show simulation results for four other
lasers, including colasing systems, with different erbium concentrations.

To show that the proposed model can indeed be used with different erbium
concentrations, it was first applied to an unpublished result obtained for a
2815 nm fiber laser pumped at 981 nm. The erbium-doped fluoride fiber had an
erbium concentration of 5 mol.%, a core radius of 5 $\mu$m, a cladding radius
of 125 $\mu$m and a NA of 0.188. The active fiber length was 14.5 m, the low
reflectivity ouput coupler FBG (LRFBG) had a reflectance of 7.1% at 2815 nm
and the high reflective input coupler FBG (HRFBG) had a reflectance of 96.5%
at the same wavelength. The output end of the laser was terminated with an
AlF3 endcap. The simulation results are showing very good agreement with
experimental data as can be seen in figure 8.

Figure 8: Simulation results of an unpublished 5 mol.% laser using the
proposed model and the WI and SI parameters. The attenuation at the laser
wavelength was fixed to 50 dB/km. The pump wavelength was 981 nm, and the
laser wavelength was 2815 nm.

In obtaining the results of figure 8, equations (24) to (26) were used without
any modification, i.e. based on the very same set of parameters.

Following the work by Guo et al. [12], we attempted to reproduce the
experimental results from the laser in ref. [2] by taking the ESA at 1.6
$\mu$m from level 1 to 3 into account. Note that this laser used a
significantly lower erbium concentration of 1 mol.%. Our simulation results
along with experimental data are shown in figure 9.

Figure 9: Simulation results of the 2825 nm signal of the laser from [2] using
the proposed model and the WI and SI parameters. The attenuation was fixed to
30 dB/km at 2825 nm. The pump wavelength was fixed to 976 nm. Loss values for
all the components were set to the values stated in ref. [2] without any
modification.

Results for the 2825 nm signal are for this case in less good agreement with
the experimental data, as the laser efficiency is a bit too low. Also, at such
low erbium concentrations, the energy transfer model has a much lower impact
on predicted laser power than for the 5 mol.%, 6 mol.% and 7 mol.% lasers,
which can be seen on figure 9. The WI model and the proposed model both give
very similar results.

A colasing system using a 1.5 mol.% fiber from ref. [8] was also simulated.
The main difference with the other studied lasers in this work is the use of
broadband mirrors instead of the narrowband FBGs. For such lasers, a redshift
of the laser wavelength with increasing pump power is usually observed [11,
9]. Our simulation results along with experimental data are shown in figure
10.

Figure 10: Simulation results of the 2.8 $\mu$m signal of the laser from [8]
using the proposed model and the WI and SI parameters. The attenuation
spectrum for FiberLabs monomode fluoride fibers was used, which is 17 dB/km at
976 nm, 30 dB/km at 1613 nm and 43 dB/km at 2887 nm. The pump wavelength was
fixed to 976 nm. Loss values for all the components were set to the values
stated in the article without any modification. The filter loss was set to
20%.

We observed good agreement between our predicted laser powers and the
experimental results. The proposed model uses a spectrum of evenly distributed
wavelength channels, so it was possible to observe a redshift in the output
spectrum of the simulation. We observed wavelengths of 2871 nm at 4 W of pump
power and 2887 nm at higher pump powers. The predicted redshift vs. pump power
is much stronger in our simulation than in reality, which could be due to
error in the measured cross section spectra of the 2.8 $\mu$m laser
transition, which in turn affects the peak gain value at each pump power.

Finally, as shown in figure 11, we also reproduced results from the laser in
ref. [7] which correspond to a 2.8 $\mu$m laser also involving broadband
mirrors and a 6 mol.% erbium doped fiber.

Figure 11: Simulation of the laser from [7] using the proposed model and the
WI and SI parameters. The attenuation was fixed to 0.55 dB/m at 2887 nm. The
pump wavelength is 975 nm, and the laser wavelength ranges from 2802 nm to
2887 nm. Intracavity signal re-injection loss of 40% and pump injection loss
of around 40% on both sides of the cavity were used.

We needed to consider very high intracavity losses to properly simulate this
laser, and no information on these losses can be found in ref. [7]. The
simulated lasing wavelength was 2802 nm at 1 W of pump power, 2824 nm at 20,
30 and 40 W, 2871 nm at 50 and 60 W and 2887 nm at higher pump powers. It
cannot be compared to ref. [7] because no information on the precise lasing
wavelength is given.

###  III-C Transfer rate analysis

Figure 4 shows that the proposed model predicts signal powers that are similar
to those predicted with the WI model at low pump powers, while it is much more
accurate at higher pump powers. Accordingly, it is interesting to follow the
ETU1 transfer rate, which has a direct impact on lasing efficiency, of the
proposed model as compared to that of the WI model. This comparison is
illustrated in figure 12 where we show how the factor in front of $N_{1}^{2}$
changes with the input pump power, whereas figure 13 shows the corresponding
atomic population evolutions. Equation (16) is used for the WI model and
equation (24) is used for the proposed model.

Figure 12: Factor in front of $N_{1}^{2}$ for ETU1 in the proposed model and
the WI model vs launched pump power at the fiber entry (pump side) calculated
using the simulation of the laser shown in figure 4. Atomic populations of
levels 0 and 1 were taken from figure 13. Figure 13: Atomic population vs.
pump power at the fiber entry (pump side) for the simulation results with the
proposed model shown in figure 4.

The factor in front of $N_{1}^{2}$ is close to the WI value at low pump power,
which is expected since the resulting signal powers in figure 4 for both
models are close at such low powers. However, the factor quickly rises with
pump power to reach a near constant value of about 2 times the WI parameter,
which explains why the proposed model leads to higher signal powers at these
pump power values. This behavior of the transfer rate is attributed to the
shape of the $N_{1}$ curve shown in figure 13.

The proposed model can also easily be used to add energy transfer processes to
the atomic population rate equations, as long as the relevant cross sections
to compute the overlap integrals in equation (22) are known. For example, the
cross relaxation processes CR3 and CR6, which are the opposite of ETU1 and
ETU2 respectively, were tested. These were tested before in ref. [11], but the
transfer rates were unknown. In our case, they could easily be computed using
the overlap integrals between $\sigma_{31}$ and $\sigma_{01}$ for CR3 and
between $\sigma_{62}$ and $\sigma_{02}$ for CR6. The McCumber relation was
used to get the cross sections that were unknown. The overlap value is
$1\times 10^{-59}$ m5 for CR3 and $1\times 10^{-57}$ m5 for CR6. However, our
simulations showed that their effect on laser efficiency is negligible since
the atomic populations of levels 3 and 6 are very low which confirms the
results from Gorjan et al. [11].

##  IV Conclusion

In summary, we have shown that the WI and SI models which were proposed to
model high-power erbium-doped fluoride fiber lasers do not reproduce well the
highly doped fiber laser powers and efficiencies. We proposed a cubic
dependence of the energy transfer rate and found that combining both Dexter’s
model and the hopping model of migration assisted energy transfer gave the
best results. Eight high-power erbium-doped fluoride fiber lasers operating
around 2.8 $\mu$m were simulated to prove that the proposed model gives better
agreement with experimental data for multiple erbium concentrations. The
efficiency drop which is generally observed at high pump powers was reproduced
with good accuracy, which means the developed model can be used for accurate
power scaling predictions. However, there is still place for improvement. For
example, we believe it should be possible to consider more general energy
transfer equations which are not the combination of the extreme cases of
energy transfer without migration and purely migration assisted energy
transfer. We believe the proposed model will help designing the next
generation of erbium-doped fluoride fiber and the cavity parameters needed to
produce average power in excess of 100W around 2.8 microns.

## Acknowledgments

The authors would like to thank Fonds de recherche du Québec – Nature et
technologies (144 616, CO256655), Canada Foundation for Innovation (5180) and
Natural Sciences and Engineering Research Council of Canada (IRCPJ469414-13,
RGPIN-2016- 05877) for financing this work, as well as Vincent Fortin and
Sébastien Magnan-Saucier for helpful discussions and technical assistance.

## References

  * [1] F. Jobin et al., “Recent developments in lanthanide-doped mid-infrared fluoride fiber lasers [Invited],” Opt. Express, vol. 30, no. 6, p. 8615-8640, Mar. 2022, doi: 10.1364/OE.450929.
  * [2] Y. O. Aydın et al., “Diode-pumped mid-infrared fiber laser with 50% slope efficiency,” Optica, vol. 4, no. 2, p. 235-238, Feb. 2017, doi: 10.1364/OPTICA.4.000235.
  * [3] D. Faucher, M. Bernier, G. Androz, N. Caron, and R. Vallée, “20 W passively cooled single-mode all-fiber laser at 2.8$\mu$m,” Opt. Lett. 36, 1104-1106, Apr. 2011, doi: 10.1364/OL.36.001104
  * [4] Y. O. Aydin, V. Fortin, R. Vallée and M. Bernier, “Towards power scaling of 2.8 $\mu$m fiber lasers,” Opt. Lett., vol. 43, no. 18, p. 4542-4545, Sep. 2018, doi: 10.1364/OL.43.004542.
  * [5] M. Bernier, D. Faucher, N. Caron, and R. Vallée, “Highly stable and efficient erbium-doped 2.8 $\mu$m all fiber laser,” Opt. Express, vol. 17, no. 19, p. 16941-16946, Sep. 2009, doi: 10.1364/OE.17.016941.
  * [6] V. Fortin, M. Bernier, S. T. Bah, and R. Vallée, “30 W fluoride glass all-fiber laser at 2.94$\mu$m” Opt. Lett., vol. 40, no. 12, p. 2882-2885, Jun. 2015, doi: 10.1364/OL.40.002882.
  * [7] S. Tokita, M. Murakami, S. Shimizu, M. Hashida, and S. Sakabe, “Liquid-cooled 24 W mid-infrared Er:ZBLAN fiber laser,” Opt. Lett., vol. 34, no. 20, p. 3062-3064, Oct. 2009, doi: 10.1364/OL.34.003062.
  * [8] J. Li, L. Wang, H. Luo, J. Xie, and Y. Liu, “High Power Cascaded Erbium Doped Fluoride Fiber Laser at Room Temperature,” IEEE Photon. Technol. Lett., vol. 28, no. 6, pp. 673–676, Mar. 2016, doi: 10.1109/LPT.2015.2504462.
  * [9] J. Li and S. D. Jackson, “Numerical Modeling and Optimization of Diode Pumped Heavily-Erbium-Doped Fluoride Fiber Lasers,” IEEE J. Quantum Electron., vol. 48, no. 4, pp. 454–464, Apr. 2012, doi: 10.1109/JQE.2012.2183856.
  * [10] M. Pollnau and S. D. Jackson, “Energy recycling versus lifetime quenching in erbium-doped 3-$\mu$m fiber lasers,” IEEE J. Quantum Electron., vol. 38, no. 2, pp. 162–169, Feb. 2002, doi: 10.1109/3.980268.
  * [11] M. Gorjan, M. Marincek, and M. Copic, “Role of Interionic Processes in the Efficiency and Operation of Erbium-Doped Fluoride Fiber Lasers,” IEEE J. Quantum Electron., vol. 47, no. 2, pp. 262–273, Feb. 2011, doi: 10.1109/JQE.2010.2073447.
  * [12] C. Guo et al., “Theoretical Analysis of ESA-Enhanced 2.8 $\mu$m Lasing in Er-Doped ZBLAN Fiber Lasers,” J. Lightwave Technol., vol. 40, no. 13, pp. 4397–4414, Jul. 2022, doi: 10.1109/JLT.2022.3163272.
  * [13] R. I. Woodward and M. Gorjan, “Modeling mid-infrared fiber laser systems,” in Mid-Infrared Fiber Photonics, Elsevier, 2022, pp. 743–801. doi: 10.1016/B978-0-12-818017-4.00003-3.
  * [14] W. J. C. Grant, “Role of Rate Equations in the Theory of Luminescent Energy Transfer,” Phys. Rev. B, vol. 4, no. 2, pp. 648–663, Jul. 1971, doi: 10.1103/PhysRevB.4.648.
  * [15] L. Agazzi, K. Wörhoff, and M. Pollnau, “Energy-Transfer-Upconversion Models, Their Applicability and Breakdown in the Presence of Spectroscopically Distinct Ion Classes: A Case Study in Amorphous Al2 O3:Er3+,” J. Phys. Chem. C, vol. 117, no. 13, pp. 6759–6776, Apr. 2013, doi: 10.1021/jp4011839.
  * [16] D. A. Zubenko, M. A. Noginov, V. A. Smirnov, and I. A. Shcherbakov, “Different mechanisms of nonlinear quenching of luminescence,” Phys. Rev. B, vol. 55, no. 14, pp. 8881–8886, Apr. 1997, doi: 10.1103/PhysRevB.55.8881.
  * [17] D. L. Dexter, “A Theory of Sensitized Luminescence in Solids,” J. Chem. Phys, vol. 21, no. 5, pp. 836–850, May 1953, doi: 10.1063/1.1699044.
  * [18] B. Di Bartolo, Ed., Energy Transfer Processes in Condensed Matter, vol. 114. in NATO ASI Series B: Physics, vol. 114. Boston, MA: Springer US, 1984, pp. 156 - 180, doi: 10.1007/978-1-4613-2407-2.
  * [19] S. Jackson, “Energy transfer processes in rare-earth-doped glass fiber,” in Mid-Infrared Fiber Photonics, Elsevier, 2022, pp. 305–331. doi: 10.1016/B978-0-12-818017-4.00004-5.
  * [20] B. Srinivasan, R. K. Jain, and G. Monnom, “Indirect measurement of the magnitude of ion clustering at high doping densities in Er:ZBLAN fibers,” J. Opt. Soc. Am. B, vol. 17, no. 2, p. 178-181, Feb. 2000, doi: 10.1364/JOSAB.17.000178.
  * [21] R. S. Quimby, W. J. Miniscalco, and B. A. Thompson, “Excited-state absorption at 980 nm in erbium-doped glass,” presented at the Fiber Laser Sources and Amplifiers III, M. J. F. Digonnet and E. Snitzer, Eds., Boston, United States, Feb. 1992, pp. 72–79. doi: 10.1117/12.134972.
  * [22] D. E. McCumber, “Einstein Relations Connecting Broadband Emission and Absorption Spectra,” Phys. Rev., vol. 136, no. 4A, pp. A954–A957, Nov. 1964, doi: 10.1103/PhysRev.136.A954.
  * [23] R. S. Quimby, “Range of validity of McCumber theory in relating absorption and emission cross sections,” Journal of Applied Physics, vol. 92, no. 1, pp. 180–187, Jul. 2002, doi: 10.1063/1.1485112.
  * [24] Y. D. Huang, M. Mortier, and F. Auzel, “Stark level analysis for Er3+-doped ZBLAN glass,” Optical Materials, vol. 17, no. 4, pp. 501–511, Sep. 2001, doi: 10.1016/S0925-3467(01)00109-4.
  * [25] W. Miniscalco, “Optical and Electronic Properties of Rare Earth Ions in Glasses,” in Rare-Earth-Doped Fiber Lasers and Amplifiers, Revised and Expanded, M. Digonnet, Ed., in Optical Science and Engineering, vol. 71. CRC Press, 2001. doi: 10.1201/9780203904657.ch2.
  * [26] B. Wang et al., “Excited state absorption cross sections of 4I13/2 of Er3+ in ZBLAN,” Optical Materials, vol. 31, no. 11, pp. 1658–1662, Sep. 2009, doi: 10.1016/j.optmat.2009.03.015.
  * [27] U. Ganiel, A. Hardy, G. Neumann, and D. Treves, “Amplified spontaneous emission and signal amplification in dye-laser systems,” IEEE J. Quantum Electron., vol. 11, no. 11, pp. 881–892, Nov. 1975, doi: 10.1109/JQE.1975.1068543.
  * [28] J. Wang and D. J. Simkin, “Energy transfer within a regular distribution of donors and acceptors: Application to the upconversion dynamics of Er3+: YAlO3 and Er3+:YAG,” Phys. Rev. B, vol. 52, no. 5, pp. 3309–3314, Aug. 1995, doi: 10.1103/PhysRevB.52.3309.
  * [29] Le Verre Fluoré, “Catalog 2023”, leverrefluore.com. Accessed: Aug. 10, 2023. [Online]. Available: https://leverrefluore.com/wp-content/uploads/2023/07/LVF-Catalog-2023.pdf

|  William Bisson received a B. Sc. degree in physics from Laval University,
Quebec City, Qc, Canada in 2022. He is currently pursuing his M. Sc. in
physics at the Center of Optics, Photonics and lasers. His research interests
include modeling of mid-infrared erbium-doped fluoride fiber lasers and
neodymium silica fiber lasers.
---|---
|  Alexandre Michaud received a B. Eng. degree in physical engineering from
Laval University, Quebec City, Qc, Canada in 2022. He is currently pursuing
his M. Sc. in physics at the Center of Optics, Photonics and lasers. His
research interests include the development of Raman fiber lasers in the mid-
infrared and neodymium silica fiber lasers.
---|---
|  Pascal Paradis received a B. Eng. degree in engineering physics, a masters
degree and a Ph. D. degree in Physics from Laval University, Quebec City, Qc,
Canada in 2017, 2019 and 2023. He is currently a postdoctoral fellow at the
Center for Optics, Photonics and Lasers, Laval University. His research
interests include pulsed mid-infrared fiber lasers for spectroscopy, remote
sensing, material processing and biomedical applications.
---|---
|  Réal Vallée received his Ph.D. from Laval University in 1986. He was a
postdoctoral fellow at the Laboratory for Laser Energetics, at the University
of Rochester before returning to Laval University as a professor in the
Department of Physics. In 2000 he was appointed Director of the Center for
Optics Photonics and Lasers (COPL), the province of Quebec network center of
excellence in photonics, and held this position until 2020. His research
interests are in fiberoptics components and their applications, non-linear
propagation of ultrafast pulses in fiber, inscription of waveguides with
femtosecond pulses and the study of infrared glasses for integrated optics. He
currently holds an NSERC-funded Industrial Research Chair in Femtosecond
Photo-Inscribed Photonic Components and Devices. Professor Vallée has
supervised over 100 graduate students, has authored over 200 peer-reviewed
papers in high impact journals and holds 16 patents. He successfully
spearheaded a grant application to the Canada Foundation for Innovation for
the construction of a building entirely dedicated to optics and photonics
university research and training. He was promoted Fellow of the Optical
Society of America in 2011 and is the recipient of the Brockhouse and Urgel-
Archambault prizes.
---|---
|  Martin Bernier received the Ph.D. degree in 2010 in physics from Laval
University, Quebec City, QC, Canada, where he is currently Full Professor with
the Department of Physics, Engineering Physics and Optics. His research
projects conducted with the Center for Optics, Photonics and Lasers, Laval
University, involve writing of fiber Bragg gratings using femtosecond pulses
and their application to the development of innovative fiber-based components,
particularly in the field of fiber lasers and sensors.
---|---