On the generality of symmetry breaking and dissipative freezing in quantum
trajectories

Joseph Tindall1$\star$, Dieter Jaksch3,2,1 and Carlos Sánchez Muñoz4

1 Clarendon Laboratory, University of Oxford, United Kingdom

2 The Hamburg Centre for Ultrafast Imaging, Hamburg, Germany

3 Institut für Laserphysik, Universität Hamburg, Hamburg, Germany

4 Departamento de Física Teórica de la Materia Condensada and Condensed Matter
Physics Center (IFIMAC), Universidad Autónoma de Madrid, 28049 Madrid, Spain

⋆<EMAIL_ADDRESS>
## Abstract

Recently, several studies involving open quantum systems which possess a
strong symmetry have observed that every individual trajectory in the Monte
Carlo unravelling of the master equation will dynamically select a specific
symmetry sector to ‘freeze’ into in the long-time limit. This phenomenon has
been termed ‘dissipative freezing’, and in this paper we demonstrate how it is
a general consequence of the presence of a strong symmetry in an open system
and independent of any microscopic details. Specifically, through several
simple, plausible, mathematical arguments, we detail how this phenomenon will
happen in most setups with only a few exceptions. Using a number of example
systems we illustrate these arguments, highlighting the relationship between
the spectral properties of the Liouvillian in off-diagonal symmetry sectors
and the time it takes for freezing to occur. In the limiting case that
eigenmodes with purely imaginary eigenvalues are manifest in these sectors,
freezing fails to occur. Such modes indicate the preservation of information
and coherences between symmetry sectors of the system and are associated with
phenomena such as non-stationarity and synchronisation. The absence of
freezing at the level of a single quantum trajectory provides a simple,
computationally efficient way of identifying these traceless modes.

###### Contents

  1. 1 Introduction
  2. 2 Theory
    1. 2.1 The GSKL Equation and the Strong Symmetry
    2. 2.2 The Monte-Carlo Unravelling of the GSKL Equation and Dissipative Freezing
    3. 2.3 Arguments for the Emergence of Dissipative Freezing
    4. 2.4 Exceptions to Dissipative Freezing
  3. 3 Results and Discussion
    1. 3.1 Numerical Details
    2. 3.2 Example: Random Matrices
    3. 3.3 Example: Coupled Qudits
    4. 3.4 Example: Lossy Bosonic Chain
  4. 4 Conclusion and Outlook

##  1 Introduction

Symmetries are fundamental to our understanding of nature. Quantum physics is
no exception to this, and the invariance of the Hamiltonian of a closed
quantum system under a given transformation can have far reaching implications
which include eigenstate degeneracy [1], an absence of thermalisation [2, 3]
and the protection of topological order [4, 5].

In an open quantum system, the influence of an external environment
complicates matters and must be accounted for when identifying symmetries and
assessing their impact on the system. The last decade has seen significant
progress in this regard, with extensive classifications of the symmetries of
open systems having been achieved [6, 7, 8, 9, 10] alongside an understanding
of their physical consequences in terms of steady state degeneracy [11, 12],
transport properties [13, 14, 15], non-stationarity [10, 16], quantum
synchronisation [17, 18], off-diagonal long-range order [19, 20] and
dissipative phase transitions [21, 22, 23]. Results such as this have made it
clear that the combination of symmetries and an external environment provides
a new pathway for realising unique, exploitable states of matter in a quantum
system.

Recently, a new symmetry-based phenomenon, termed ‘dissipative freezing’, was
witnessed in a large quantum spin precessing under the influence of
dissipation parallel to the axis of rotation [24]. It was observed how the
individual quantum trajectories - solutions to the stochastic Monte Carlo
unravelling of the master equation - randomly selected a specific symmetry
sector to converge into in the long-time limit; breaking that symmetry at the
trajectory level. Due to the simplicity of the model at hand an analytical
proof of this unexpected result was provided, although it was unclear whether
this phenomenon could be observed in more general setups. Since this result,
several independent works have noted the manifestation of this freezing in
other open system setups [25, 26], but without offering a deeper explanation
for its origin.

In this work, we unify these results and demonstrate that dissipative freezing
is a general feature of open quantum systems with a strong symmetry. We
achieve this by providing plausible mathematical arguments for its occurrence,
alongside identifying the connection between this phenomenon and previous
mathematical works establishing the ergodicity of open systems when the
external environment is subject to continuous monitoring. We then proceed to
give several illustrative numerical examples, from a range of open systems,
where dissipative freezing can be observed — relating the results of our
simulations to our earlier mathematical ones.

Importantly, the understanding of dissipative freezing provided in this work
also allows us identify the exceptional situations in which, despite the
presence of a strong symmetry, dissipative freezing is absent. The most
prominent of these situations involves the closing of the spectral gap in the
off-diagonal symmetry subspaces and heralds the emergence of traceless long-
lived modes in the Liouvillian’s spectrum. These modes ensure the preservation
of coherences and information between symmetry sectors despite the
environmental influence and are associated with complex phenomena such as non-
stationarity and synchronisation [10, 17, 18]. By framing this gap closing in
terms of an absense of dissipative freezing we provide a novel,
computationally efficient method for identifying the presence of such
traceless modes in open systems.

##  2 Theory

###  2.1 The GSKL Equation and the Strong Symmetry

Our starting point is a quantum system coupled to an environment under the
Markov approximation. The most general completely positive trace-preserving
(CPTP) map for the evolution of this system’s density matrix $\rho(t)$ at time
$t$ is given by the Gorini-Sudarshan-Kosakowski-Lindblad (GSKL) equation [27,
28]:

$\frac{\partial\rho(t)}{\partial t}=\mathcal{L}\rho(t)=-{\rm
i}[H,\rho(t)]+\sum_{j=1}^{M}\gamma_{j}\left[L_{j}\rho(t)L_{j}^{\dagger}-\frac{1}{2}\\{L^{\dagger}_{j}L_{j},\rho(t)\\}\right],$
(1)

where $H$ is the system’s Hamiltonian and $L_{1},L_{2},...,L_{M}$ are a series
of ‘jump’ operators which describe the interaction of the system with the
environment. We set $\hbar=1$ throughout and use $\gamma_{j}$ to set the
dissipation strength.

Due to the addition of the dissipative term in the equation of motion of the
system symmetries can no longer simply be identified by operators which
commute with the Hamiltonian $H$ and there are, in fact, several different
types of symmetry that an open system can possess [6]. Here, we focus on the
‘strong symmetry’, which can be identified by a Hermitian operator $A$ that
satisfies

$[H,A]=[L_{j},A]=[L_{j}^{\dagger},A]=0\ \forall j,$ (2)

ensuring the observable $\langle A\rangle$ is a constant of motion.

As $A$ is Hermitian, we can diagonalise it

$A=\sum_{\alpha=1}^{D}\lambda_{\alpha}\sum_{\beta=1}^{d_{\alpha}}\ket{\alpha,\beta}\bra{\alpha,\beta},$
(3)

where the index $\alpha$ runs over the $D$ distinct eigenvalues
$\lambda_{\alpha}$ of $A$ and $\beta$ runs over the corresponding $d_{\alpha}$
degenerate eigenvectors $\ket{\alpha,\beta}$ for a given $\alpha$, i.e.
$A\ket{\alpha,\beta}=\lambda_{\alpha}\ket{\alpha,\beta}$. Throughout this work
we will use $\alpha$ to notate a ‘subspace’ or ‘symmetry subspace’ spanned by
these degenerate vectors, and define a ‘block’ as the operator space spanned
by the $d_{\alpha}d_{\alpha^{\prime}}$ elements
$\ket{\alpha,\beta}\bra{\alpha^{\prime},\beta^{\prime}}$ with
$\beta=1...d_{\alpha}$ and $\beta^{\prime}=1...d_{\alpha^{\prime}}$. We refer
to a given block as ‘diagonal’ if $\alpha=\alpha^{\prime}$ and ‘off-diagonal’
otherwise.

Now, as the Hamiltonian and the jump operators all commute with $A$ it
immediately follows that
$A(O\ket{\alpha,\beta})=\lambda_{\alpha}(O\ket{\alpha,\beta})$ for
$O=H,L_{j},L_{j}^{\dagger}\ \forall j$ and so they cannot map an eigenstate of
$A$ out of its subspace. They are thus all ‘block-diagonal’, meaning we can
write them in the following form

$O=\sum_{\alpha=1}^{D}\sum_{\beta,\beta^{\prime}=1}^{d_{\alpha}}O_{\alpha,\beta,\beta^{\prime}}\ket{\alpha,\beta}\bra{\alpha,\beta^{\prime}}.$
(4)

The presence of a strong symmetry therefore immediately implies there exists a
simultaneous block-decomposition of the Hamiltonian and jump operators. We
emphasize that the converse is also true, if there exists some basis
$\\{\ket{\alpha,\beta}\\}$ in which the Hamiltonian and jump operators are
simultaneously block-diagonal then there exists a strong symmetry of the
system given by the operator $A$ in Eq. (3) — i.e. the operator which is just
some scalar multiple of the identity matrix in each diagonal block and zero
everywhere else.

Due to this strong symmetry there will be multiple steady states of the
system, i.e. multiple density matrices $\rho_{\infty}$ which satisfy
$\mathcal{L}\rho_{\infty}=0$ [6]. In fact, there will be at least $D$ trace
unity Hermitian steady states. This can be understood from the fact that the
projection of the Liouvillian into the superoperator space formed from two
copies of a given diagonal block constitutes a valid CPTP map and that via
Evan’s theorem there will exist at least one trace unity steady state for such
a map [29, 30]. There can also exist traceless steady states, or even
traceless imaginary eigenmodes (modes $\rho_{\rm imag}$ with imaginary
eigenvalues — i.e. $\mathcal{L}\rho_{\rm imag}={\rm i}\lambda\rho_{\rm imag},\
\lambda\in\mathbb{R}_{\neq 0}$) of the Liouvillian which originate from the
projection of the Liouvillian into the superoperator subspace formed from two
copies of a given off-diagonal block. These eigenmodes involve coherences
between states in different symmetry sectors, i.e. they have ${\rm
Tr}(\rho\ket{\alpha,\beta}\bra{\alpha^{\prime},\beta^{\prime}})\neq 0$ for
$\alpha\neq\alpha^{\prime}$. We will refer to them collectively as traceless,
non-decaying modes and later in this work we will demonstrate how such modes
prevent dissipative freezing.

###  2.2 The Monte-Carlo Unravelling of the GSKL Equation and Dissipative
Freezing

In order to fully introduce dissipative freezing we must describe the Monte-
Carlo unravelling of the GSKL equation in Eq. (1). This is an unravelling of
the master equation in Eq. (1) into a stochastic equation of motion at the
level of pure states [31, 32, 33]. After statistical averaging of the
independent solutions, known as quantum trajectories, to this equation of
motion the exact dynamics of the density matrix of the system is guaranteed to
be recovered. We emphasize that this unravelling of the GSKL equation is not
just a mathematical formalism: the stochastic equation of motion which governs
the trajectories exactly describes the dynamics of the open system when the
environment is subject to continuous measurement with an appropriate choice of
POVM (positive-operator-valued-measure) [34, 35]. As a result, the phenomenon
of dissipative freezing that we will describe is relevant from both a
mathematical and a physical perspective.

Figure 1: Dissipative Freezing: The Hilbert space of an open quantum system
fragments into a series of subspaces in the presence of a strong symmetry
operator $A$. For an initial state in a superposition of states in these
different subspaces the long-time dynamics of a single quantum trajectory will
almost always involve a ‘freezing’ into one of these subspaces selected at
random.

Given an initial pure state of the system, $\ket{\psi(0)}$, all independent
solutions to the stochastic equation of motion which unravels Eq. (1) can be
written, to first order in the propagation timestep $\Delta t$, in the
following form

$\ket{\psi^{(i)}(t)}=\frac{1}{\mathcal{N}_{i}}\prod_{k=1}^{n}X^{(i)}_{k}\ket{\psi(0)}$
(5)

where $i$ is used to index the given trajectory, $\mathcal{N}_{i}$ is a
normalisation constant, $t=n\Delta t$ is the time with
$n\in\mathbb{Z}^{+}_{0}$ and

$X^{(i)}_{k}\in S=\left\\{\sqrt{\gamma_{1}\Delta
t}L_{1},...,\sqrt{\gamma_{M}\Delta t}L_{M},1-{\rm i}H_{\rm eff}\Delta
t\right\\}=\\{S_{1},...,S_{M+1}\\},$ (6)

with $H_{\rm eff}=H-\frac{\rm i}{2}\sum_{j}\gamma_{j}L^{\dagger}_{j}L_{j}$.
The state $\ket{\psi^{(i)}(n\Delta t)}$ is proportional to
$X^{(i)}_{n}\ket{\psi^{(i)}((n-1)\Delta t)}$ and the operator $X^{(i)}_{n}$
which propagated it from time $(n-1)\Delta t$ to time $n\Delta t$ is
determined stochastically from the set $S$ with the respective ‘draw’
probabilities being

$p_{m}=\begin{cases}\langle\psi^{(i)}((n-1)\Delta
t)|S^{\dagger}_{m}S_{m}|\psi^{(i)}((n-1)\Delta t)\rangle\qquad&m\leq M,\\\
1-\sum_{m=1}^{M}\langle\psi^{(i)}((n-1)\Delta
t)|S^{\dagger}_{m}S_{m}|\psi^{(i)}((n-1)\Delta
t)\rangle\qquad&m=M+1.\end{cases}$ (7)

The evolution of a single trajectory is thus a Markovian process: the choice
of operator $X^{(i)}_{n}$ is dependent only on the previous state
$\ket{\psi^{(i)}((n-1)\Delta t)}$.

The density matrix of the system at time $t$ — i.e. the solution to Eq. (1) —
can be approximated as a sum over $N$ independent trajectories:
$\rho(t)\approx\frac{1}{N}\sum_{i=1}^{N}\ket{\psi^{(i)}(t)}\bra{\psi^{(i)}(t)}$
and in the limits $\Delta t\rightarrow 0$ and $N\rightarrow\infty$ this
approximation becomes exact [32]. The unravelling of Eq. (1) we have just
described is accurate to first order in $\Delta t$. We opt for this first
order unravelling in this paper as it is illustrative and all of our arguments
are valid for arbitrarily small $\Delta t$. Unravellings which are to higher
order in $\Delta t$ involve letting the wavefunction evolve freely under ${\rm
exp}(-{\rm i}H_{\rm eff}\Delta t)$ until its norm decays below a certain
threshold, at which point it is acted on with a randomly chosen jump operator.

Now that we have introduced the Monte Carlo unravelling of the GSKL equation
we can explicitly describe the phenomenon of dissipative freezing:

Dissipative Freezing: If an open system governed by the GSKL equation
possesses a strong symmetry $A$ then, regardless of the initial state of the
system, in the limit $t\rightarrow\infty$ a given quantum trajectory will,
generically (i.e. outside of just a few exceptions), only have non-zero
coefficients in one symmetry subspace.

We should emphasize that the given subspace into which freezing occurs will
vary with each trajectory, allowing the system to recover the full density
matrix — which could have non-zero coefficients in multiple diagonal blocks —
upon averaging. Figure 1 provides an illustration of this phenomenon.

In the case when the initial state of the system is already confined to one
symmetry subspace then the phenomenon of dissipative freezing is immediate and
trivial. Each trajectory will always remain in that subspace as the operators
$X^{(i)}_{k}$ are block-diagonal and cannot map it out of the given subspace.
This also means that once a trajectory is frozen, it will remain frozen for
any future times. When the initial state of the system is in a superposition
of states in different subspaces then dissipative freezing is non-trivial. It
implies a dynamical breaking of the strong symmetry at the level of
trajectories — with the system randomly selecting one of the (possibly)
degenerate subspaces of $A$ in the long-time limit. As we will see, such a
symmetry breaking is a direct consequence of the presence of a strong symmetry
and means that the system preserves the symmetry at the ensemble level by
breaking it at the level of pure states. In the following section, we will
provide several mathematical arguments for the occurrence of this phenomenon.
These arguments will also allow us to identify the exceptional cases where
freezing does not occur.

###  2.3 Arguments for the Emergence of Dissipative Freezing

#### Recovering a Block-diagonal Steady State

The simplest argument for dissipative freezing can be understood from the form
of the steady state. Assuming there are no traceless non-decaying eigenmodes
of the Liouvillian, then this steady state can always be written in block
diagonal form

$\rho_{\infty}=\lim_{t\rightarrow\infty}\exp(\mathcal{L}t)\rho(0)=\sum_{\alpha=1}^{D}\sum_{\beta,\beta^{\prime}=1}^{d_{\alpha}}\rho_{\alpha,\beta,\beta^{\prime}}\ket{\alpha,\beta}\bra{\alpha,\beta^{\prime}},$
(8)

just like the jump operators and Hamiltonian. We know that, if $\Delta
t\rightarrow 0$ in our trajectories unravelling, upon averaging over an
infinite number of quantum trajectories in the long-time limit we must exactly
recover this steady state, i.e.

$\rho_{\infty}=\lim_{N\rightarrow\infty}\lim_{t\rightarrow\infty}\frac{1}{N}\sum_{i=1}^{N}\ket{\psi^{(i)}(t)}\bra{\psi^{(i)}(t)},\qquad\ket{\psi^{(i)}(t)}=\sum_{\alpha,\beta}c^{(i)}_{\alpha,\beta}(t)\ket{\alpha,\beta},$
(9)

where the $c^{(i)}_{\alpha,\beta}(t)$ are the coefficients of the $i$th
trajectory at time $t$ for the basis vector $\ket{\alpha,\beta}$. Equations
(8) and (9) therefore combine to put a constraint on the coefficients of the
trajectories in the long-time limit via

$\lim_{N\rightarrow\infty}\lim_{t\rightarrow\infty}\frac{1}{N}\frac{\partial^{m}}{\partial
t^{m}}\sum_{i=1}^{N}c^{(i)}_{\alpha,\beta}(t)c^{*\
(i)}_{\alpha^{\prime},\beta^{\prime}}(t)=0,\qquad\forall\alpha\neq\alpha^{\prime}\qquad
m\in\mathbb{Z}^{+}_{0},$ (10)

ensuring the steady state is, and remains, block-diagonal in the long-time
limit. It is clear that if dissipative freezing occurs then Eq. (10) is
satisfied as, for each $i$, at least one of the coefficients in the product
$c^{(i)}_{\alpha,\beta}(t)c^{*\ (i)}_{\alpha^{\prime},\beta^{\prime}}(t)$
(with $\alpha\neq\alpha^{\prime}$) will always be $0$ and remain $0$ for
$t\rightarrow\infty$; ensuring the product, and all its time-derivatives,
vanish for each $i$.

Without dissipative freezing, Eq. (10) requires the product
$c^{(i)}_{\alpha,\beta}(t)c^{*\ (i)}_{\alpha^{\prime},\beta^{\prime}}(t)$ —
and all its time-derivatives — to ensemble average111In this paper we define
ensemble-average as the stochastic averaging over all trajectories. to $0$
without the individual terms being $0$. The product of trajectory coefficients
in the same subspaces must also ensemble average to the steady state elements,
i.e. $(1/N)\sum_{i}c^{(i)}_{\alpha,\beta}(t)c^{*\
(i)}_{\alpha,\beta^{\prime}}(t)=\rho_{\alpha,\beta,\beta^{\prime}}$ as
$t\rightarrow\infty$. Imposing Eq. (10) on top of this becomes very
restrictive without dissipative freezing, placing constraints on the product
of coefficients in different subspaces — which can have a very different
structure — alongside the standard constraints on coefficients in the same
subspace. Moreover, from an abstract perspective, we can continuously add
additional symmetry subspaces to the system (by performing a direct sum of the
Hamiltonian and jump operators with some additional operators) — each with
their own associated steady state — and, without dissipative freezing, Eq.
(10) would impose ever-increasing constraints on the coefficients. Dissipative
freezing thus provides a natural mechanism by which this restriction is
satisfied independently of the specific details and number of symmetry
subspaces.

#### Long-Products of Matrices

For our second argument, and an illustration of how dissipative freezing
emerges as a dynamical process, we consider the explicit time evolution of a
single trajectory. To help with our analysis, we introduce the projection
operator
$P_{\alpha}=\sum_{\beta=1}^{d_{\alpha}}\ket{\alpha,\beta}\bra{\alpha,\beta}$
which projects into the subspace spanned by the eigenvectors of the strong
symmetry $A$ with eigenvalue $\lambda_{\alpha}$. For a given trajectory at
time $t$ we can form the projection into the $\alpha$ subspace as
$\ket{\psi_{\alpha}^{(i)}(t)}=P_{\alpha}\ket{\psi^{(i)}(t)}$. Due to the
block-diagonal nature of the Hamiltonian and jump operators it follows from
Eq. (5) that

$\ket{\psi_{\alpha}^{(i)}(t)}=\frac{1}{\mathcal{N}_{i}}\prod_{k=1}^{n}X^{(i)}_{\alpha,k}\ket{\psi_{\alpha}(0)},$
(11)

where $X^{(i)}_{\alpha,k}=P_{\alpha}X^{(i)}_{k}P_{\alpha}$. The norm of the
wavefunction in Eq. (11) can then be interpreted as the probability
$p^{(i)}(\alpha,t)$ of the trajectory $i$ being observed in the subspace
$\alpha$ at time $t$. We scale this by the normalisation constant
$\mathcal{N}_{i}$ (which is the same for each subspace and so irrelevant to
our analysis) to define the un-normalised weight

$w^{(i)}(\alpha,t)=\bra{\psi_{\alpha}(0)}(A^{(i)}_{\alpha}(t))^{\dagger}A^{(i)}_{\alpha}(t)\ket{\psi_{\alpha}(0)},\qquad
A^{(i)}_{\alpha}(t)=\prod_{k=1}^{n}X^{(i)}_{\alpha,k}.$ (12)

In the limit $t\rightarrow\infty$ Eq. (12) is the expectation value of an
infinitely long product of matrices. Long products of matrices have been
studied extensively in mathematics and physics due to their relevance to
dynamical systems and chaos [36, 37, 38]. When building up repeated matrix
products such as
$A^{(i)}_{\alpha}(t)=X^{(i)}_{\alpha,n}X^{(i)}_{\alpha,n-1}X^{(i)}_{\alpha,n-2}...X^{(i)}_{\alpha,1}$,
the quantities of interest are usually the singular values
$\sigma_{\alpha,q}^{(i)}(n)\coloneqq\sigma_{q}(A^{(i)}_{\alpha}(t))$ of the
product with $q$ indexing them from largest to smallest. These singular values
are intrinsically related to the powers of the matrix norm
$||A^{(i)}_{\alpha}(t)||^{q}$ whose logarithm is often known as the $q$th
Lyaponuv exponent. In general, such quantities values will grow or decay
exponentially with $n$ [39, 40], causing the largest singular value
$\sigma^{(i)}_{\alpha,1}(n)$ to become dominant when taking expectation values
and allowing us to write

$(A_{\alpha}^{(i)}(t))^{\dagger}A^{(i)}_{\alpha}(t)\approx\big{(}\sigma^{(i)}_{\alpha,1}(n)\big{)}^{2}P(\sigma^{(i)}_{\alpha,1}(n)),$
(13)

where $P(\sigma^{(i)}_{\alpha,1}(n))$ is a projector to the subspace (not to
be confused with the subspaces indexed by $\alpha$) spanned by the (possibly
degenerate) eigenvectors of
$(A^{(i)}_{\alpha}(t))^{\dagger}A^{(i)}_{\alpha}(t)$ with the associated
eigenvalue $\big{(}\sigma^{(i)}_{\alpha,1}(n)\big{)}^{2}$. Provided the
initial state in a given pair of symmetry subspaces $\alpha$ and
$\alpha^{\prime}$ has a non-zero weight in their respective highest singular-
value subspaces then the ratio of the probabilities $p^{(i)}(\alpha,t)$ and
$p^{(i)}(\alpha^{\prime},t)$ will be proportional to the ratios of the
corresponding singular values, i.e.

$\frac{w^{(i)}(\alpha,t)}{w^{(i)}(\alpha^{\prime},t)}=\frac{p^{(i)}(\alpha,t)}{p^{(i)}(\alpha^{\prime},t)}\propto\left(\frac{\sigma^{(i)}_{\alpha,1}(n)}{\sigma^{(i)}_{\alpha^{\prime},1}(n)}\right)^{2}.$
(14)

If these singular values each vary exponentially with $n$, however, then the
ratio of probabilities will clearly either vanish or diverge in the limit
$t\rightarrow\infty$ depending on whether the growth rate of the leading
singular value is largest in $\alpha$ or $\alpha^{\prime}$. One symmetry
subspace thus becomes completely dominant and the wavefunction effectively
vanishes in the other subspaces. Freezing thus occurs between all symmetry
subspaces with distinct growth rates as a natural consequence of the
exponentially growing nature of $A^{(i)}_{\alpha}(t)$.

We should emphasize that our statement of an exponential growth in the leading
singular values $\sigma^{(i)}_{\alpha,q}(n)$ is intended to be that of an
average change by orders of magnitude in linear time. We have been careful not
to assign a specific growth rate to this $\forall n$ and instead talk about
the ‘average’ change as we are working with matrices and only under certain
conditions can one assign a single growth rate $\forall n$ to the leading
singular value of long products [36]. Nonetheless, the average change of
orders of magnitude in linear time will, via Eq. (14), be reflected in the
$p^{(i)}(\alpha,t)$ and cause them to differ by increasing orders of magnitude
as time increases linearly. In our numerical results we will observe this
growth explicitly.

#### Ergodicity and Quantum Systems Subject to Repeated Measurement

For a final perspective behind dissipative freezing, we will discuss its
connection between established mathematical results on the ergodicity of open
systems. Specifically, work in the early 2000s successfully proved the
ergodicity of an open quantum system undergoing continuous measurement in the
environment when there is a single unique steady state [41, 42]. This
ergodicity guarantees that, as $t\rightarrow\infty$ the time-average of a
single trajectory can be used to recover the ensemble average of multiple
trajectories and thus the dynamics of the open system when left to evolve
freely - without continuous monitoring.

In this work it was also realised that, if the steady state of the system is
dependent on the initial state in some manner, then this ergodic theorem
cannot be proven [41]. Such a result connects to dissipative freezing as the
presence of a strong symmetry guarantees the existence of multiple steady
states and thus a dependence of
$\lim_{t\rightarrow\infty}\exp(\mathcal{L}t)\rho(0)$ on the initial state
$\rho(0)$ via its weights in the different subspaces. An ergodic theorem thus
cannot be established in such a case — i.e. the time-averaged dynamics of a
single trajectory cannot be used to recover the ensemble average over multiple
trajectories. This is clearly connected to our earlier discussion on the
restrictiveness of Eq. (10).

This lack of a standard ergodic theorem for a system with degenerate steady
states implies that something else must be at play in order for the ensemble
average over trajectories to recover the dynamics of the open system. This was
uncovered in Ref. [43], where Kuemmerer and Maassen proved that, in this
degenerate case, a single trajectory will establish ergodicity solely with
respect to one of the, randomly selected, steady states. By observing that the
existence of a strong symmetry in an open system implies multiple steady
states and a series of symmetry subspaces the work of Ref. [43] provides an
ergodic understanding of the phenomenon of dissipative freezing. By freezing
into one of the symmetry subspaces, the trajectory of an open system can
establish ergodicity and — upon time averaging — recover the steady state
solely within that subspace. If the probability of freezing into a specific
subspace is proportional to the weight of the initial state in that subspace
then the ensemble average over trajectories will recover the steady state
specific to the given initial state and ensure the conservation of $\langle
A\rangle$ at the ensemble level.

###  2.4 Exceptions to Dissipative Freezing

There are ‘exceptional’ cases where the conservation of $\langle A\rangle$ and
the recovery of the appropriate steady state is achieved within the
trajectories formalism without dissipative freezing occurring. Here we
identify the two cases we are aware of and which follow naturally from the
mathematical arguments of Section 2.3.

#### Similar Symmetry Subspaces

If there exist two subspaces $\alpha_{1}$ and $\alpha_{2}$ where the following
is true

$P_{\alpha_{1}}HP_{\alpha_{1}}=P_{\alpha_{2}}HP_{\alpha_{2}},\qquad
P_{\alpha_{1}}L_{j}P_{\alpha_{1}}=\exp{({\rm
i}\theta_{j})}P_{\alpha_{2}}L_{j}P_{\alpha_{2}},\ \forall
j\qquad\theta_{j}\in\mathbb{R},$ (15)

then freezing cannot occur between them. This is straightforwardly proven by
observing that
$(A^{(i)}_{\alpha_{1}}(t))^{\dagger}A^{(i)}_{\alpha_{1}}(t)=(A^{(i)}_{\alpha_{2}}(t))^{\dagger}A^{(i)}_{\alpha_{2}}(t)$
at all times and thus the weights in the two subspaces cannot diverge with
respect to one another. We term the symmetry subspaces ‘similar’ in this case.

#### Traceless Modes

The second, less trivial, case occurs when there are traceless non-decaying
eigenstates of the Liouvillian. Such states possess coherences
$\ket{\alpha,\beta}\bra{\alpha^{\prime},\beta^{\prime}}$ between states in
different symmetry subspaces, $\alpha\neq\alpha^{\prime}$ and do not decay
away in the long-time limit. These non-block-diagonal modes directly imply
that Eq. (10) cannot be true (it hinges on the assumption that
$\lim_{t\rightarrow\infty}\exp(\mathcal{L}t)\rho(0)$ is always block-diagonal)
and as dissipative freezing always ensures the validity of this equation it is
clear these traceless modes imply the absence of dissipative freezing.
Physically, these traceless modes respresent information and coherences
between symmetry subspaces which is protected from the system’s interaction
with the environment. These modes play an important role in the emergence of
non-stationarity in open systems, [10], quantum synchronisation [17, 18] and
quantum information processing [44].

In the following section we will illustrate dissipative freezing in a number
of numerical examples of open quantum systems. We will also observe the
exceptions we have just described.

##  3 Results and Discussion

###  3.1 Numerical Details

For all of the following examples we evolve our open system using the first-
order trajectories routine described in Section 2.2. To minimise any incurred
errors in the trajectory routine we set $\Delta t\omega=0.001$ throughout,
where $\omega$ is the central energy scale in the Hamiltonian and will be
specified for each example. We always initialise the system in an equal
superposition of states in the symmetry sectors of ‘interest’ — i.e. the
weight of the initial state is the same for each symmetry sector where it is
non-zero. In our numerics, we also calculate the ‘freeze-time’ of a given
trajectory as the time $t$ it takes for the probability $p(\alpha^{*},t)$ to
exceed a threshold value $1-\epsilon$ in one of the subspaces $\alpha^{*}$. If
freezing is going to occur, $p(\alpha^{*},t)$ is the probability, at time $t$,
the trajectory will freeze to the subspace $\alpha^{*}$ as
$t\rightarrow\infty$. Here we set $\epsilon=10^{-10}$.

###  3.2 Example: Random Matrices

As our first example we consider the Lindblad equation, as in Eq. (1), with
Hamiltonian and Lindblad operators which are block-diagonal random matrices —
i.e. the off-diagonal blocks are set to $0$ and the diagonal blocks of both
the Hamiltonian and Lindblad operators consist of matrices drawn from the
space of random Hermitian matrices and the complex Ginibre ensemble [45]
respectively. This open system is then guaranteed to possess a strong symmetry
in the form of an operator which is a different scalar multiple of the
identity in each block.

In Fig. 2a-b we plot the evolution of the weights (both normalised and un-
normalised) of a single trajectory in each symmetry subspace, along with the
corresponding leading singular value of the matrix $A_{\alpha}(t)$ defined in
Eq. (12) — dropping the $i$ superscript here as we are focussed solely on a
single trajectory. We use a logarithmic $y$-axis in order to expose the
exponential growth of these quantities. Freezing occurs into the subspace
which has the highest overall growth rate and the normalised weight in this
subspace saturates to unity in the long-time limit. The exponential growth in
the unnormalised weights coincides with an exponential growth of the leading
singular values of the evolution matrix $A_{\alpha}(t)$ in each subspace
$\alpha$. We also observe that the ratio of the largest and second largest
singular values associated with this matrix is increasing at an exponential
rate justifying the approximation in Eq. (13). This growth is intrinsically
tied to the initial state being propagated under an increasingly long-product
of matrices.

Figure 2: Freezing for a ‘random’ Liouvillian constructed from a random block-
diagonal Hamiltonian and a random block-diagonal jump operator. There are four
blocks, each of dimension $4$ and, in each block, the Hamiltonian $H$ consists
of a random Hermitian matrix scaled by a factor of $\omega$ whilst the single
jump operator consists of a random operator drawn from the complex Ginibre
ensemble [45]. We set $\gamma=4\omega$. The system is initialised in a
normalised superposition of random states in each subspace with the initial
weight in each subspace the same. a) Logarithm of the unormalised weight of
the wavefunction in each subspace versus time for a single trajectory. Inset)
Normalised weight of the wavefunction versus time. b) Logarithm of the largest
singular values of the matrix $A_{\alpha}(t)$ in each subspace $\alpha$ for
the same single trajectory. Inset) Logarithm of the ratio of the largest and
second largest singular value in each subspace — values are not plotted once
they can no longer be accurately determined. c) Probability density function
for the freeze-time $t\omega$ for $10^{4}$ trajectories. The two different
colours correspond to two independent realisations of the random Liouvillian —
the initial state remains the same between them. Inset) Cumulative density
function for the freeze time.

In Fig. 2c we then plot histograms of the freeze times when measured over
$10^{4}$ trajectories for two different realisations of the random Liouvillian
— where a given realisation is constructed by drawing the block-diagonal
Hamiltonian and jump operator from their respective ensembles. The two
histograms of freeze-times display qualitatively different distributions even
though the matrices are drawn from the same underlying ensemble both times.
This is due to the fact that the freeze-time will be strongly controlled by
the part of the eigenspectrum corresponding to the ‘off-diagonal’ modes of the
Liouvillian, i.e. the eigenstates of the Liouvillian which are not block-
diagonal and contain coherences between states in different symmetry
subspaces. Such coherences must manifest in the transient dynamics of the
system and so the system cannot freeze until all these off-diagonal modes have
died away sufficiently. The spacing between the eigenvalues of the longest-
lived of these modes will thus play a crucial factor in the freezing
statistics.

In this example this spacing will clearly change for each draw of the block-
matrices from their ensembles. In the next example we will consider a non-
random, more physical Liouvillian, allowing us to more straightforwardly
calculate the eigenvalues for these off-diagonal, low-lying modes and directly
compare them with the corresponding freeze-times.

###  3.3 Example: Coupled Qudits

For this example we move onto a more physical setup consisting of two coupled
qudits — $a$ and $b$ — with qudit $a$ undergoing continuous magnetic dephasing
in the $z$ direction. The Hamiltonian reads

$H=\omega(s^{x}_{a}s^{x}_{b}+s^{y}_{a}s^{y}_{b}+s^{z}_{a}s^{z}_{b}),$ (16)

where $s^{\alpha}_{a/b}$ is the spin-3/2 operator acting in the $\alpha=x,y,z$
direction on either qubit $a$ or $b$. The dephasing is described by the single
jump operator $L=s^{z}_{a}$.

This system possesses a strong symmetry in the form of the total $z$
magnetisation $A=s^{z}_{a}+s^{z}_{b}$ and thus the Hamiltonian and jump-
operator can be simultaneously block-diagonalised into seven subspaces of
dimensions $1,2,3,4,3,2,1$ and with respective total $z$ magnetisations
$-3/2,-1,-1/2,0,1/2,1$ and $3/2$. Whilst there are seven subspaces, we should
not expect freezing between all of them. Specifically, for each pair of
subspaces with the same absolute magnetisation Eq. (15) is valid and they are
sufficiently similar to prevent freezing. We emphasize that there are no non-
decaying traceless modes of the Liouvillian here and Eq. (10) holds for these
pair of subspaces even though dissipative freezing does not occur. This can be
understood from the similarity of the subspaces: apart from the initial state,
the only difference in the dynamics in the two subspaces is that when a jump
occurs the coefficients of the trajectory in the negative magnetisation
subspace pick up a factor of $e^{{\rm i}\pi}$ compared to the positive one.
After a certain time the total phase picked up will thus be $n\pi$, where $n$
is the number of jumps. The summation in the expression
$\sum_{i}c^{(i)}_{\alpha^{\prime},\beta^{\prime}}(t)c^{(i)}_{\alpha^{\prime\prime},\beta^{\prime\prime}}(t)$
can thus be split into two, one involving trajectories where $n$ is odd and
one where $n$ is even. In the long-time limit the probability that $n$ is even
or odd should be equal and so the two terms should be identical aside from a
scale factor of $e^{{\rm i}\pi}$, which causes them to cancel out and the
ensemble-average — and all it’s derivatives — to equal to $0$ without freezing
occurring.

Figure 3: Freezing for coupled qudits governed by the GSKL equation with the
Hamiltonian defined in Eq. (16) and a single jump operator $L=s^{z}_{a}$. We
set $\gamma=3\omega$. The index $\alpha=1,2,..7$ refers to subspaces of total
$z$ magnetisation $-3/2,-1,...,1,3/2$ respectively. a) Logarithm of the un-
normalised weight of a single trajectory — initialised in a superposition of
random states in each subspace — versus time in each of the magnetisation
subspaces. Inset) Normalised weight for $t\omega\leq 50$. b) Probability
density of the freeze-time $t$ for $10^{4}$ trajectories and for trajectories
initialised in a superposition of states in the indicated subspaces with
$\gamma=3\omega$. c) Average freeze-time (black dots) vs $\gamma\omega$ on a
log-log scale for trajectories initialised in a superposition of random states
in the $\alpha=2$ and $\alpha=3$ subspaces respectively. Rosebrown curve shows
the function $-28.1/\Delta$ where $\Delta$ is the inter-sector spectral gap
for these subspaces — this gap has been calculated by exact diagonalisation of
the Liouvillian for the full range of values of $\gamma\omega$.

In Fig. 3 we plot the evolution of trajectories in this system and illustrate
these results, observing freezing between subspaces of different absolute
magnetisation due to the distinct exponential growth of the matrices
$A_{\alpha}(t)$ whilst observing an absence of freezing in subspaces with
magnetisation of opposite sign due to the similarity of the subspaces. We then
calculate the freezing times for trajectories initialised in a superposition
of two different symmetry sectors. We observe distinctive distributions for
each symmetry sector, indicating that they are specific to the form of the
matrices in those sectors. We select the $\alpha=2$ and $\alpha=3$ sectors and
analyse this further, varying $\gamma$ and computing the average freezing time
and the inter-sector spectral gap $\Delta$. This ‘inter-sector spectral gap’
is calculated as the smallest real part of the eigenvalues of the Liouvillian
whose corresponding eigenmodes contain coherences between states in the
relevant symmetry sectors (e.g. here these sectors are $\alpha=2$ and
$\alpha=3$).

For small values of $\gamma$ we find both the inverse spectral gap and the
freezing-time scale linearly with the inverse of $\gamma$. They are thus
proportional for small $\gamma\omega$, with the proportionality constant shown
in Fig. 3 obtained by numerical fitting and dependent on more complex features
of the Liouvillian’s spectrum. It is clear from our results that as the gap
closes the freezing time should diverge and thus freezing will be absent. This
should be completely general as the closing of the spectral gap means that Eq.
(10) cannot hold — the Liouvillian possesses traceless modes with eigenvalues
having real part $0$ and so the general form of the steady state cannot be
block-diagonal. In the next example we will observe this explicitly with a
system whose gap closes for finite $\gamma$.

For intermediate values of $\gamma$ the relationship between the freeze-time
and the spectral gap is not so clear, with identical values of the spectral
gap resulting in substantially different freeze-times. This is likely due to
the effect of other eigenvalues lying close to the spectral gap on the
freezing statistics, with the freeze-time being dependent on this distribution
and not just the spectral gap. In the small $\gamma$ regime this low-lying
distribution will simply scale with $\gamma$, allowing us to ‘ignore’ the
distribution’s effect on the freeze-times and focus solely on it’s distance to
the imaginary axis. For systems where the lowest-lying inter-sector mode is
well separated from the rest of the spectrum, we expect that the inter-sector
gap can be related to the freezing time beyond just the small $\gamma$ regime.

###  3.4 Example: Lossy Bosonic Chain

For our final example we consider a periodic boundary chain of $L$ sites
hosting non-interacting bosons collectively coupled to a cavity. In momentum
space, the Hamiltonian reads

$H=\omega a^{\dagger}a-g(a+a^{\dagger})\sum_{k}b_{k}^{\dagger}b_{(k+\pi)\ {\rm
mod}\ 2\pi}-2J\sum_{k}\cos(k)n_{k},$ (17)

where $a$ and $a^{\dagger}$, respectively annihilate and create photons in the
cavity whilst $b_{k}$ and $b^{\dagger}_{k}$ respectively annihilate and create
bosons in the chain with momenta
$k\in\\{\frac{2\pi}{L},\frac{4\pi}{L},...,\frac{2\pi L}{L}\\}$. We also define
$n_{k}$ as the boson number operator for momentum $k$, i.e.
$n_{k}=b_{k}^{\dagger}b_{k}$. Photon loss is introduced into the cavity via a
single jump operator $L=a$ and we fix ourselves to half-filling, i.e. the
total number of bosons is always $L/2$. The master equation formed from these
operators was studied in Ref. [26] in the context of dissipative phase
transitions of a cavity-immersed quantum gas. Importantly this setup possesses
a series of $L/2$ commuting strong symmetry operators of the form
$S_{k}=n_{k}+n_{k+\pi}$ with $k\leq\pi$.

In Ref. [26] it was observed that, for the initial states considered, the
individual trajectories in this system always freeze to a symmetry subspace.
Here we unravel this phenomenon of dissipative freezing using the intuition we
have built up. Notably we also describe how, beyond a certain system size and
for certain initial states, dissipative freezing is not guaranteed between all
subspaces due to the presence of traceless non-decaying modes in the
Liouvillian.

First, as the strong symmetry operators all commute there clearly exists a
simultaneous block-decomposition of the Hamiltonian and jump-operator.
Specifically, the corresponding subspaces or blocks can be indexed by the
unique tuples $(s_{1},s_{2},...,s_{L/2})$ where $s_{i}$ corresponds to the
number of bosons in momentum modes $k=\frac{2\pi i}{L}$ and $k=\frac{2\pi
i}{L}+\pi$, i.e. $s_{i}=\langle n_{k}+n_{k+\pi}\rangle$. If the total number
of bosons is $N$ then it follows the number of unique tuples (or subspaces) is
$D=\frac{(N+L/2-1)!}{N!(L/2-1)!}$ [46] which grows exponentially with $L$ for
$N=L/2$.

An initial state initialised in a superposition of states with different
tuples is therefore expected to undergo the phenomenon of dissipative freezing
— which is exactly what we observe in Fig. 4a-b for a small $L=4$ site chain
with the initial state having a non-zero weight in all subspaces. The
characteristic exponential growth in the weights and singular values is
observed just like in our earlier examples. Moreover, for small values of
$\gamma\omega$ we can again find the freezing times are directly proportional
to the inverse of the inter-sector spectral gap.

It might be natural to assume that freezing will occur in this setup for all
system sizes. In Fig. 4c-d we treat a system for $L=6$ and observe that this
is not the case, with freezing only occurring between certain subspaces. The
origin of this lack of freezing is traceless modes of the Liouvillian, with
eigenvalue $0$, which appear for $L=6$. In our numerics we also observe that
these modes also appear for $L>6$ and independently of the maximum number of
photons that we enforce. For $L=6$ these modes connect the symmetry sectors
where $s_{1}+s_{2}$ is the same, which likely is related to the fact that the
$\cos(k)$ term in the Hamiltonian only differs by a minus sign for the two
separate (i.e. not simply related by a shift of $+\pi$) momentum modes
$k=2\pi/3$ and $k^{\prime}=4\pi/3$ which correspond to $s_{1}$ and $s_{2}$
respectively. When $L\geq 6$ we can always identify such pairs of modes where
$\cos(2\pi i/L)=-\cos(2\pi j/L)$ and $i,j\in{1,2,...,L/2}$. For $L=4$ this is
not possible.

Figure 4: Freezing dynamics for non-interacting bosons collectively coupled to
a single-mode electromagnetic field. The Hamiltonian is defined in Eq. (17)
and there is a single jump operator $L=\sqrt{\gamma}a$. We set
$\gamma=5\omega$, $g=2\omega$, $J=2\omega$ and the system is always
initialised in an equal superposition of random states in each symmetry
subspace. The index $\alpha$ runs over the valid momentum tuples
$(s_{1},s_{2},...,s_{L/2})$ which specify the symmetry subspaces, the explicit
relationship between $\alpha$ and these tuples is detailed in appendix A. We
enforce a maximum photon number in our simulations of $n_{max}=5$. In plots
a-c) we set $L=4$ whilst in plots d-e) we set $L=6$. a) Logarithm of the
weight of a single trajectory in each of the strong symmetry subspaces verus
time. b) Logarithm of the largest singular value of $A_{\alpha}(t)$ in each
subspace versus time. c) Average freeze-out time on a log-log plot for
trajectories initialised in a superposition of states in subspaces $\alpha=1$
and $\alpha=3$. Rosebrown curve shows $-11.0/\Delta$ where $\Delta$ is the
inter-sector spectral gap for $\alpha=2$ and $\alpha=3$ — which has been
calculated by exact diagonalisation of the Liouvillian for the full range of
values of $\gamma\omega$ . d) Logarithm of the weight of a single trajectory
in each of the strong symmetry subspaces versus time. Inset) Dynamics of the
normalised weight for $t\omega\leq 40$. e) Absolute values of the elements of
the matrix $p(\alpha,t)p(\alpha^{\prime},t)$ for $t\omega=200$ and averaged
over $10^{4}$ trajectories.

These traceless modes are clearly manifest as non-zero off-diagonal elements
in the ensemble-averaged matrix $p(\alpha,t)p(\alpha^{\prime},t)$ plotted in
Fig. 4e. As we know that there are no similar subspaces (i.e. pairs of
subspaces where Eq. (15) holds) in our setup we can be confident this matrix
provides a clear picture of where the traceless non-decaying modes of the
system reside. Even a single trajectory initialised across all subspaces,
however, allows us to infer this information, because if there are traceless
modes and the freezing time diverges then any given trajectory should have a
negligible probability of freezing on any finite time-scale.

Figure 4d explicitly demonstrates this information being inferred at the
single-trajectory level: all pairs of subspaces which share a traceless mode
have weights of the same order of magnitude for all time, in direct agreement
with the ensemble averaged picture provided by Fig. 4e. This single-trajectory
picture thus immediately points to a steady state which is not block-diagonal
for initial states containing coherences which overlap with these modes. The
freezing time thus diverges but, as opposed to the previous setups, this
happens for all finite $\gamma\omega$ and does not require it to to tend to
$0$ or $\infty$.

Importantly, our results in this paper show we do not need to wait until
$t\rightarrow\infty$ to infer this lack of freezing at the single trajectory
level. Figure 4 pictures a clear manifestation of these traceless modes after
a finite period of time, with the probabilities in the associated subspaces
saturating to a constant, non-zero non-unity value. If the trajectory is left
to evolve over a time-scale sufficiently longer than that associated with any
decay modes of the Liouvillian, then any signatures of freezing (e.g.
divergent weights in different subspaces) should be apparent. If they are not
then it is reasonable to assume freezing will not occur and traceless non-
decaying modes are present. Given the identification of a strong symmetry and
the corresponding symmetry subspaces, the ability to utilise a single
trajectory to detect the existence of these traceless modes is much more
computationally efficient method than having to resort to ensemble averaging
or memory-intensive calculations in the superoperator picture.

We note that these traceless modes and the resultant absence of freezing
between certain subspaces was not observed in Ref. [26] as the initial states
used did not span the necessary subspaces to have an overlap with these modes.

##  4 Conclusion and Outlook

In this work we have demonstrated that the recently-observed phenomenon of
dissipative freezing is a general consequence of the presence of a strong
symmetry in an open system. The associated symmetry-breaking in the
trajectories picture is the way the system ensures symmetry is preserved at
the ensemble level. We have also introduced a range of examples where
dissipative freezing can be directly observed and quantified.

In these examples, we showed how the freezing statistics are dependent on the
lowest lying eigenvalue in the ‘off-diagonal’ superoperator subspaces. When
the real part of this eigenvalue vanishes then the Liouvillian posseses
traceless non-decaying modes, immediately implying an absence of dissipative
freezing for any single given trajectory. We have thus identified a
computationally efficient method for identifying the existence of such modes.

We envisage a number of future avenues of research stemming from this work.
Firstly, identifying whether further exceptions to this phenomenon are
possible or whether those identified in this work — traceless non-decaying
modes and similar subspaces — are the only possible exceptions would greatly
further our understanding of the effect of symmetries on trajectories in open
quantum systems.

Secondly, rigorously quantifying how the spectral distribution of the lowest-
lying inter-sector Liouvillian eigenvalues affects the freezing statistics of
the system could lead to novel ways to determine said portion of the spectrum
of open systems from within the memory efficient trajectories formalism.

Lastly, and by no means least, an important further pursuit would be to
quantify the role that symmetry-breaking perturbations play in the breakdown
of dissipative freezing. Such perturbations would lead to an opening of the
Liouvillian’s gap in the full superoperator space and quantifying this through
the freezing statistics could provide new insights into dissipative phase
transitions in open systems.

##### Acknowledgements

C.S.M. acknowledges that the project which gave rise to these results received
the support of a fellowship from "la Caixa" foundation (ID 100010434) and from
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie Grant Agreement No. 847648, with fellowship code
LCF/BQ/P120/11760026, and financial support from the Proyecto Sinérgico CAM
2020 Y2020/TCS-6545 (NanoQuCo-CM). DJ and JT acknowledge funding from EPSRC
grant EP/P009565/1. DJ also acknowledges funding from the Cluster of
Excellence ‘Advanced Imaging of Matter’ of the Deutsche Forschungsgemeinschaft
(DFG) - EXC 2056 - project ID 390715994. We are also grateful to Berislav Buča
for discussions on dissipative freezing and for suggesting the coupled qudit
model.

## Appendix A: Indexing the subspaces of the lossy bosonic chain

In Table 1 we tabulate the relationship between the single subspace index
$\alpha$ we adopt in the figures and the unique tuples
$(s_{1},s_{2},...,s_{L/2})$ which explicitly define the strong symmetry
subspaces for the lossy bosonic chain in Section 3.4.

$\mathbf{L=4}$:  $(s_{1},s_{2})$ $\alpha$ (0,2) 1 (1,1) 2 (2,0) 3 |  $\mathbf{L=6}$:  $(s_{1},s_{2},s_{3})$ $\alpha$ (0,0, 3) 1 (0,1,2) 2 (0,2, 1) 3 (0,3,0) 4 (1,0,2) 5 (1,1,1) 6 (1,2,0) 7 (2,0,1) 8 (2,1,0) 9 (3,0,0) 10
---|---

Table 1: Relationship between the subspace index $\alpha$ we adopt in Figure 4
and the momentum tuples ($s_{1},s_{2},....,s_{L/2}$) which define the symmetry
subspaces in the lossy bosonic chain. Left table is for $L=4$ and right is for
$L=6$.

## References

  * [1] D. J. Gross,  _The role of symmetry in fundamental physics_ ,  Proceedings of the National Academy of Sciences 93(25), 14256 (1996),  10.1073/pnas.93.25.14256,  https://www.pnas.org/content/93/25/14256.full.pdf.
  * [2] Y. Tang, W. Kao, K.-Y. Li, S. Seo, K. Mallayya, M. Rigol, S. Gopalakrishnan and B. L. Lev,  _Thermalization near integrability in a dipolar quantum Newton’s cradle_ ,  Phys. Rev. X 8, 021030 (2018),  10.1103/PhysRevX.8.021030.
  * [3] M. Rigol and M. Fitzpatrick,  _Initial-state dependence of the quench dynamics in integrable quantum systems_ ,  Phys. Rev. A 84, 033640 (2011),  10.1103/PhysRevA.84.033640.
  * [4] Z.-C. Gu and X.-G. Wen,  _Tensor-entanglement-filtering renormalization approach and symmetry-protected topological order_ ,  Phys. Rev. B 80, 155131 (2009),  10.1103/PhysRevB.80.155131.
  * [5] Y. You, T. Devakul, F. J. Burnell and S. L. Sondhi,  _Subsystem symmetry protected topological order_ ,  Phys. Rev. B 98, 035112 (2018),  10.1103/PhysRevB.98.035112.
  * [6] B. Buča and T. Prosen,  _A note on symmetry reductions of the Lindblad equation: transport in constrained open spin chains_ ,  New Journal of Physics 14(7), 073007 (2012),  10.1088/1367-2630/14/7/073007.
  * [7] V. V. Albert and L. Jiang,  _Symmetries and conserved quantities in Lindblad master equations_ ,  Phys. Rev. A 89, 022118 (2014),  10.1103/PhysRevA.89.022118.
  * [8] S. Lieu, M. McGinley and N. R. Cooper,  _Tenfold way for quadratic Lindbladians_ ,  Phys. Rev. Lett. 124, 040401 (2020),  10.1103/PhysRevLett.124.040401.
  * [9] A. Altland, M. Fleischhauer and S. Diehl,  _Symmetry classes of open fermionic quantum matter_ ,  Phys. Rev. X 11, 021037 (2021),  10.1103/PhysRevX.11.021037.
  * [10] B. Buča, J. Tindall and D. Jaksch,  _Non-stationary coherent quantum many-body dynamics through dissipation_ ,  Nature Communications 10(1), 1730 (2019),  10.1038/s41467-019-09757-y.
  * [11] Z. Zhang, J. Tindall, J. Mur-Petit, D. Jaksch and B. Buča,  _Stationary state degeneracy of open quantum systems with non-abelian symmetries_ ,  Journal of Physics A: Mathematical and Theoretical 53(21), 215304 (2020),  10.1088/1751-8121/ab88e3.
  * [12] P. Zhang and Y. Chen,  _Violation and revival of Kramers’ degeneracy in open quantum systems_ (2021), 2108.05493.
  * [13] D. Manzano and P. I. Hurtado,  _Symmetry and the thermodynamics of currents in open quantum systems_ ,  Phys. Rev. B 90, 125138 (2014),  10.1103/PhysRevB.90.125138.
  * [14] V. Popkov and R. Livi,  _Manipulating energy and spin currents in non-equilibrium systems of interacting qubits_ ,  New Journal of Physics 15(2), 023030 (2013),  10.1088/1367-2630/15/2/023030.
  * [15] J. Thingna, D. Manzano and J. Cao,  _Dynamical signatures of molecular symmetries in nonequilibrium quantum transport_ ,  Scientific Reports 6(1), 28027 (2016),  10.1038/srep28027.
  * [16] F. Iemini, A. Russomanno, J. Keeling, M. Schirò, M. Dalmonte and R. Fazio,  _Boundary time crystals_ ,  Phys. Rev. Lett. 121, 035301 (2018),  10.1103/PhysRevLett.121.035301.
  * [17] J. Tindall, C. S. Muñoz, B. Buča and D. Jaksch,  _Quantum synchronisation enabled by dynamical symmetries and dissipation_ ,  New Journal of Physics 22(1), 013026 (2020),  10.1088/1367-2630/ab60f5.
  * [18] B. B. Buča, C. Booker and D. Jaksch,  _Algebraic theory of quantum synchronization and limit cycles under dissipation_ (2021), 2103.01808.
  * [19] J. Tindall, B. Buča, J. R. Coulthard and D. Jaksch,  _Heating-induced long-range $\eta$ pairing in the Hubbard model_,  Phys. Rev. Lett. 123, 030603 (2019),  10.1103/PhysRevLett.123.030603.
  * [20] S. Dutta and N. R. Cooper,  _Long-range coherence and multiple steady states in a lossy qubit array_ ,  Phys. Rev. Lett. 125, 240404 (2020),  10.1103/PhysRevLett.125.240404.
  * [21] F. Minganti, A. Biella, N. Bartolo and C. Ciuti,  _Spectral theory of Liouvillians for dissipative phase transitions_ ,  Phys. Rev. A 98, 042118 (2018),  10.1103/PhysRevA.98.042118.
  * [22] M. Marcuzzi, E. Levi, W. Li, J. P. Garrahan, B. Olmos and I. Lesanovsky,  _Non-equilibrium universality in the dynamics of dissipative cold atomic gases_ ,  New Journal of Physics 17(7), 072003 (2015),  10.1088/1367-2630/17/7/072003.
  * [23] E. M. Kessler, G. Giedke, A. Imamoglu, S. F. Yelin, M. D. Lukin and J. I. Cirac,  _Dissipative phase transition in a central spin system_ ,  Phys. Rev. A 86, 012116 (2012),  10.1103/PhysRevA.86.012116.
  * [24] C. Sánchez Muñoz, B. Buča, J. Tindall, A. González-Tudela, D. Jaksch and D. Porras,  _Symmetries and conservation laws in quantum trajectories: Dissipative freezing_ ,  Phys. Rev. A 100, 042113 (2019),  10.1103/PhysRevA.100.042113.
  * [25] D. Manzano, M. A. Martínez-García and P. I. Hurtado,  _Coupled activity-current fluctuations in open quantum systems under strong symmetries_ ,  New Journal of Physics 23(7), 073044 (2021),  10.1088/1367-2630/ac0f19.
  * [26] C.-M. Halati, A. Sheikhan and C. Kollath,  _Breaking strong symmetries in dissipative quantum systems: Bosonic atoms coupled to a cavity_ ,  Phys. Rev. Research 4, L012015 (2022),  10.1103/PhysRevResearch.4.L012015.
  * [27] H. Breuer, P. Breuer, F. Petruccione and S. Petruccione,  _The Theory of Open Quantum Systems_ ,  Oxford University Press,  ISBN 9780198520634 (2002).
  * [28] Á. Rivas and S. Huelga,  _Open Quantum Systems: An Introduction_ ,  SpringerBriefs in Physics. Springer Berlin Heidelberg,  ISBN 9783642233548 (2011).
  * [29] D. E. Evans,  _Irreducible quantum dynamical semigroups_ ,  Communications in Mathematical Physics 54(3), 293 (1977),  10.1007/BF01614091.
  * [30] D. E. Evans and H. Hanche-Olsen,  _The generators of positive semigroups_ ,  Journal of Functional Analysis 32(2), 207 (1979),  https://doi.org/10.1016/0022-1236(79)90054-5.
  * [31] K. Mølmer, Y. Castin and J. Dalibard,  _Monte Carlo wave-function method in quantum optics_ ,  J. Opt. Soc. Am. B 10(3), 524 (1993),  10.1364/JOSAB.10.000524.
  * [32] A. J. Daley,  _Quantum trajectories and open many-body quantum systems_ ,  Advances in Physics 63(2), 77 (2014),  10.1080/00018732.2014.933502,  https://doi.org/10.1080/00018732.2014.933502.
  * [33] C. Gardiner and P. Zoller,  _The Quantum World of Ultra-cold Atoms and Light_ ,  Cold atoms. Imperial College Press,  ISBN 9781783264612 (2014).
  * [34] H. M. Wiseman,  _Quantum trajectories and quantum measurement theory_ ,  Quantum and Semiclassical Optics: Journal of the European Optical Society Part B 8(1), 205 (1996),  10.1088/1355-5111/8/1/015.
  * [35] T. A. Brun,  _Continuous measurements, quantum trajectories, and decoherent histories_ ,  Phys. Rev. A 61, 042107 (2000),  10.1103/PhysRevA.61.042107.
  * [36] H. Furstenberg,  _Noncommuting random products_ ,  Transactions of the American Mathematical Society 108(3), 377 (1963).
  * [37] V. I. Oseledets,  _A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems._ ,  Trans. Moscow Math. Soc. pp. 179–210 (1968).
  * [38] C. M. Newman,  _The distribution of Lyapunov exponents: Exact results for random matrices_ ,  Communications in Mathematical Physics 103(1), 121 (1986),  10.1007/BF01464284.
  * [39] J. Vanneste,  _Estimating generalized Lyapunov exponents for products of random matrices_ ,  Phys. Rev. E 81, 036701 (2010),  10.1103/PhysRevE.81.036701.
  * [40] S. Oliveira and D. E. Stewart,  _Exponential splittings of products of matrices and accurately computing singular values of long products_ ,  Linear Algebra and its Applications 309(1), 175 (2000),  https://doi.org/10.1016/S0024-3795(99)00273-6.
  * [41] J. D. Cresser,  _Ergodicity of quantum trajectory detection records_ ,  In H. J. Carmichael, R. J. Glauber and M. O. Scully, eds., _Directions in Quantum Optics_ , pp. 358–369. Springer Berlin Heidelberg, Berlin, Heidelberg,  ISBN 978-3-540-40894-9 (2001).
  * [42] B. K. mmerer and H. Maassen,  _An ergodic theorem for quantum counting processes_ ,  Journal of Physics A: Mathematical and General 36(8), 2155 (2003),  10.1088/0305-4470/36/8/312.
  * [43] B. Kümmerer and H. Maassen,  _A pathwise ergodic theorem for quantum trajectories_ ,  Journal of Physics A: Mathematical and General 37(49), 11889 (2004),  10.1088/0305-4470/37/49/008.
  * [44] D. A. Lidar, I. L. Chuang and K. B. Whaley,  _Decoherence-free subspaces for quantum computation_ ,  Phys. Rev. Lett. 81, 2594 (1998),  10.1103/PhysRevLett.81.2594.
  * [45] J. Ginibre,  _Statistical ensembles of complex, quaternion, and real matrices_ ,  J. Math. Phys. 6, 440 (1965),  10.1063/1.1704292,  https://doi.org/10.1063/1.1704292.
  * [46] P. Flajolet and R. Sedgewick,  _Analytical Combinatorics_ ,  Cambridge University Press, USA, 1 edn.,  ISBN 0521898064 (2009).


# Instances of small size with no weakly stable matching for three-sided
problem with complete cyclic preferences

Lerner E.Yu

###### Abstract

Given $n$ men, $n$ women, and $n$ dogs, each man has a complete preference
list of women, while each woman does a complete preference list of dogs, and
each dog does a complete preference list of men. We understand a matching as a
collection of $n$ nonintersecting triples, each of which contains a man, a
woman, and a dog. A matching is said to be nonstable, if one can find a man, a
woman, and a dog which belong to different triples and prefer each other to
their current partners in the corresponding triples. Otherwise the matching is
said to be stable (a weakly stable matching in 3DSM-CYC problem). According to
the conjecture proposed by Eriksson, Söstrand, and Strimling (2006), the 3DSM-
CYC problem always has a stable matching. However, Lam and Paxton have
proposed an algorithm for constructing preference lists in the 3DSM-CYC
problem of size $n=90$, which has allowed them to disprove the mentioned
conjecture. The question on the existence of counterexamples of a lesser size
remains open. In this paper, we construct an instance of the 3DSM-CYC problem
with no stable matching, whose size $n=18$. The rest counterexamples
constructed by us have the size $n=24$.

##  1 Introduction

This paper has issued from our attempts to reformulate in terms of the graph
theory the C.-K. Lam and C.G. Plaxton results [10], which disprove the
conjecture proposed by K. Eriksson, J. Söstrand, and P. Strimling [5]. Such a
formulation would make the mentioned results more demonstrative. Our paper is
the second one (after [11]) in the series of works devoted to this topic. In
the framework of the considered problem, we can treat the paper [11] as a
preparatory one; it solves a separate problem stated by P. Biró and E.
McDermid [1]. Below we give all the necessary definitions and introductory
remarks; we also prove the result from [11] that we use here. We tried to make
this paper understandable even for an unprepared reader.

Assume that there are $n$ men and $n$ women, and each one among them has a
preference list of representatives of the opposite sex. We understand a stable
matching as a partition of the set of men and women into heterosexual
families, where there exists no tuple consisting of a man and a woman, who
would prefer each other to their partners (if any). Initially, complete
preference lists were studied by D. Gale and L.S. Shaple; in [6], they prove
that a stable matching necessarily exists and propose an $O(n^{2})$-hard
algorithm for forming it. Note that the mentioned algorithm is also applicable
for the case of incomplete preference lists, but in such a case, some men and
women, possibly, remain single (single people of opposite sex do not want to
marry each other and cannot take partners from married persons). A certain
modification of the Gale–Shaple algorithm (see, e.g., [9]) allows one to find,
if possible, a matching without single men and women or to prove its absence,
otherwise.

In [9], D.E. Knuth states the question whether it is possible to generalize
the theory of stable matchings to the case of three genders. The most
interesting variant in the $k$-gender case occurs when preferences are cyclic,
i.e., representatives of the 1st gender rank representatives of the 2nd one,
the latter do representatives of the 3rd gender, etc., and each representative
of the $k$th gender has a preference list of representatives of the 1st one.

A tuple containing exactly one representative of each gender is called a
family, and a set of nonintersecting families is called a matching. A matching
is said to be nonstable, if one can find one representative of each gender
would become “happier” when forming a new family. Otherwise the matching is
said to be a weakly stable. In what follows, for brevity, we use the term “a
stable matching” instead of the term “a weakly stable matching”. Let the
number of representatives of each gender equal $n$. The authors of [4] prove
that in the case of complete preference lists, a stable matching always
exists, provided that $n\leqslant k$ (here $k$ is the number of genders). In
[5], K. Eriksson et al. generalize this result for the case when $k=3$ and
$n=k+1=4$. Ibid, they state the conjecture that the problem of finding a
stable matching in the 3-gender case with complete preference lists (three-
dimensional stable matching with cyclic preferences problem — 3DSM-CYC or just
3DSM) has a solution with any $n$. Using the statement of the satisfiability
problem and performing an extensive computer-assisted search, the authors of
[13] prove the validity of the conjecture stated by K. Eriksson et al. for
$n=5$. In [14], B. Pittel proves that with random preference lists the mean
value of stable matchings in the 3DSM problem grows as
$\Omega(n^{2}\ln^{2}(n))$.

The paper [1] by P. Biró and E. McDermid is devoted to studying the 3DSMI
problem (the 3-gender stable matching problem with incomplete preference
lists). Note that in this problem, similarly to the two-gender case, some
agents, possibly, remain single. But in contrast to the two-gender case,
according to results obtained in [1], the 3DSMI problem is not necessarily
solvable; P. Biró and E. McDermid give an explicit example of the 3DSMI
problem of size $n=6$, where no stable matching exists. Moreover, they prove
that the problem of establishing the solvability of the 3DSMI problem is NP-
complete.

The conjecture stated by K. Eriksson et al. has been recently disproved by
C.-K. Lam and C.G. Plaxton [10]. They associate the 3DSMI problem with a
certain 3DSM problem, where $n$ is 15 times greater than the initial size;
this problem is solvable if and only if so is the initial 3DSMI problem.
Therefore, the problem of establishing the solvability of the 3DSM problem is
NP-complete. Note that C.-K. Lam and C.G. Plaxton initially considered the
case of arbitrary $k>2$ (the factor 15 is a specification of the factor
$k((k-1)^{2}+1))$ for $k=3$).

The example proposed in the paper [1] allows one to construct an instance of
the 3DSM problem with no stable matching for $n=90=6\times 15$. No
counterexample of a smaller size is known by now.

An evident way to reduce the size of counterexamples of the 3DSM problem is to
solve the P. Biró and E. McDermid problem that implies the search of instances
of the 3DSMI problem with no stable matching for $n<6$. We solve the mentioned
problem in the paper [11]. We prove the absence of such instances for $n<3$
and construct several counterexamples for the 3DSMI problem with $n=3$.
Therefore, the result obtained in [10] allows one to construct an instance of
the 3DSM problem with no stable matching for $n=45$.

Another approach implies the reduction of the multiplier (its current value
equals 15) when constructing an unsolvable instance of the 3DSM problem from
an unsolvable 3DSMI problem. In this paper, we prove that one can associate
each instance of the 3DSMI problem with no stable matching of size $n$ with an
instance of the 3DSM problem with the same property of size $8n$. This allows
one to reduce the size of the counterexample of the 3DSM problem to $3\times
8=24$.

Finally, one can use constructions, which associate unsolvable 3DSMI problems
with unsolvable 3DSM problems only in certain particular cases. Developing
constructions mentioned in the previous paragraph and making use of specific
features of a certain concrete instance of the 3DSMI problem with no stable
matching, we have succeeded in proposing a counterexample for the 3DSM problem
of size $n=18$.

The rest part of the paper has the following structure. In Sect. 2, we present
formal statements of 3DSM-CYC and 3DSMI-CYC problems in terms of the graph
theory. In Sect. 3, we prove that one can associate each instance of the 3DSMI
problem of size $n$ with no stable matching with an instance of the 3DSM
problem of size $8n$ with the same property. In Sect. 4, we propose a
counterexample of the 3DSM problem of size 18. In Sect. 5, we summarize the
obtained results and discuss some potential future work. In Appendix, we give
a short proof of theorem 2 from [11] that we use in this paper.

##  2 The statement of the 3DSM (3DSMI) problem in terms of the graph theory

Let $G$ be some directed graph. Denote the set of its edges by $E(G)$; assume
that no edge is multiple. Assume that the vertex set $V(G)=V$ of the graph $G$
is divided into three subsets, namely, the set of men $M$, women $F$, and dogs
$D$. Assume that the edges $(v,v^{\prime})$, $v,v^{\prime}\in V$, of this
graph are such that either $v\in M,v^{\prime}\in F$, or $v\in F,v^{\prime}\in
D$, or $v\in D,v^{\prime}\in M$. Assume that $|M|=|F|=|D|$ (otherwise we
supplement the corresponding subgraph with vertices that are not connected
with the rest part of the graph). The number $n=|M|=|F|=|D|$ is called the
problem size. Evidently, the length of all cycles in the graph $G$ is a
multiple of $3$. Note also that this condition ensures the possibility of
dividing the vertex set of any diraph $G$ into 3 subsets $M$, $F$, and $D$ so
that all its edges are directed as is indicated above.

Each edge $(v,v^{\prime})$, $v,v^{\prime}\in V$, corresponds to some positive
integer $r(v,v^{\prime})$; it is called the rank of this edge. For fixed $v\in
V$, all possible ranks $r(v,v_{1}),\ldots,r(v,v_{k})$ coincide with
$\\{1,\ldots,k\\}$, where $k$ is the outgoing vertex degree $v$ (if
$r(v,v^{\prime})=1$, then $v^{\prime}$ is the best preference for $v$, and so
on).

We understand a three-sided matching as a subgraph $H$ of the graph $G$,
$V(H)=V(G)=V$, where each vertex $v\in V$ has at most one outgoing edge and
the following condition is fulfilled: if a vertex $v$ has an outgoing edge,
then this edge belongs to a cycle of length 3 in the graph $H$. Cycles of
length 3 in the graph $H$ are called families. Evidently, each family,
accurate to a cyclic shift, takes the form $(m,f,d)$, where $m\in M$, $f\in
F$, and $d\in D$. Note that in what follows, for convenience of denotations of
families, we do not fix the order of genders in a family, i.e., we treat
denotations of families as triples derived from an initial one by a cyclic
shift as equivalent.

In what follows, we sometimes use the notion of a family in a wider sense,
namely, as any cycle of length 3 in the graph $G$. However, if some three-
sided matching $H$ is fixed, then we describe other cycles of length 3
explicitly, applying the term “a family” only to cycles that enter in a three-
sided matching.

A matching $\mu$ is a collection of all families of a three-sided matching
$H$. For a vertex $v$, $v\in V$, in the matching $\mu$, the rank $R_{\mu}(v)$
is defined as the rank of the edge that goes out of this vertex in the
subgraph $H$. If some vertex $v$ in the subgraph $H$ has no outgoing edge,
then $R_{\mu}(v)$ is set to $+\infty$.

A triple $(v,v^{\prime},v^{\prime\prime})$ is said to be blocking for some
matching $\mu$, if it is a cycle in the graph $G$, and

$r(v,v^{\prime})<R_{\mu}(v),\quad
r(v^{\prime},v^{\prime\prime})<R_{\mu}(v^{\prime}),\quad
r(v^{\prime\prime},v)<R_{\mu}(v^{\prime\prime}).$

A matching $\mu$ is said to be stable, if no blocking triple exists for it.

The 3DSMI problem consists in finding a stable matching for a given graph $G$.
As is well known, such a matching does not necessarily exists. Moreover, the
problem of establishing its existence for a given graph $G$ is NP-complete. As
was mentioned in the introduction, the proof of this fact belongs to P. Biró
and E. McDermid. They have constructed an explicit example of the graph $G$ of
size 6, for which no stable matching exists.

The 3DSM problem represents a particular case of the 3DSMI problem, where the
outgoing (and incoming) degree of each vertex of the corresponding graph
equals the problem size $n$.

Let us construct on the set of agents (graph vertices) the same map as that
used in [12] and in other papers. Denote it by the symbol $\mu$; the same
symbol stands for the matching, from which we construct this map. If an agent
$x$ remains single, then we put $\mu(x)=x$. Otherwise $\mu(x)=y$, where $y$ is
the vertex, which represents the endpoint of an edge in the subraph $H$ that
generates the matching $\mu$. Informally speaking, in this case, $\mu(x)$ is
the partner of $x$ in the family, whom $x$ is “partial” to. Evidently, in the
3DSM problem, as distinct from the 3DSMI one, the equality $\mu(x)=x$ is
impossible for a stable matching.

In what follows, we use symbols $G$ and $H$ for graphs of problems 3DSM and
3DSMI, as well as their subgraphs; we also use, when appropriate, the symbol
$H$ with various subscripts. At the same time, when considering concurrently
the subgraph $H$ and subgraphs denoted by the symbol $H$ with various
subscripts, we distinguish them. In such cases, the symbol $H$ denotes the
“central”, in a sense, subgraph, while the symbol $H$ with some subscript does
a “peripheral” subgraph characterized by its subscript. We use the same
subscript in denotations for vertices of the latter subgraph.

##  3 The correspondence between unsolvable 3DSMI and 3DSM problems

In this section, we prove that each 3DSMI problem of size $n$ with no stable
matching corresponds to a 3DSM problem of size $8n$ with the same property. As
was mentioned in the introduction, an analogous result with the multiplier 15
belongs to C.-K. Lam and C.G. Plaxton [10]. As a corollary, making use of
results obtained by us in the paper [11] (see Lemma 3 in the Appendix), we
obtain concrete instances of the 3DSM problem of size 24 with no stable
matching.

Let us consider subgraphs of the weighted graph of the 3DSM problem, which
include some vertices and edges of the initial graph. Let us first prove the
Key Lemma 1. We “attach” one copy of the subgraph considered in this lemma to
each vertex of the graph $H$ that defines the 3DSMI problem with no stable
matching.

$c$$d$$e$$a$$t$$x$$b$$s$$r^{\prime}_{x}$$r_{x}^{\prime}+1$33 Fig. 1: The part
of the graph defining preferences that is considered in Lemma 1. Vertices of
various colors correspond to various genders. Bold lines represent edges of
rank 1, while dashed ones do those of rank 2. Ranks of edges represented by
dotted lines equal $r^{\prime}_{x}$, $r^{\prime}_{x}+1$, or 3 (ranks are
indicated near edges). Red lines represent edges that correspond to pairs that
enter in families in the matching $\mu$ in case I.2 in the proof of Lemma 1.

######  Lemma 1 (The Key Lemma)

Let some subgraph of the graph of the 3DSM problem take the form shown in Fig.
1, in particular,

$r(b,s)=r(d,t)=3,\quad r(x,c)=r^{\prime}_{x},\ r(x,d)=r^{\prime}_{x}+1,$ (1)

where $r^{\prime}_{x}\in\\{1,\ldots,n-1\\}$. Let $\mu$ be a stable matching in
this problem, while $y=\mu(x)$. Then only one of the following two
alternatives is possible:
A) $y\in\\{c,d\\}$ and $\mu(y)\in\\{a,b,t\\}$;
B) $r(x,y)<r^{\prime}_{x}$.

Proof: When proving the lemma, we repeatedly use the following fact: a triple
$(\alpha,\beta,\gamma)$, which does not belong to the mentioned matching, with
the least possible (under considered constraints) values of $r(\alpha,\beta)$,
$r(\beta,\gamma)$, and $r(\gamma,\alpha)$ is not blocking, only if one of
pairs $(\alpha,\beta)$, $(\beta,\gamma)$, or $(\gamma,\alpha)$ belongs to some
family in the matching.

I) Let us first prove the validity of Lemma 1 with $y\neq c$. Assume the
contrary. Then alternatives $r(x,y)=r^{\prime}_{x}$ and
$r(x,y)<r^{\prime}_{x}$ (see item B) are excluded. Therefore,
$r(x,y)>r^{\prime}_{x}$.

Consider the triple $(c,a,x)$. It is not blocking only in the following two
cases:
1) either $\mu(c)=a$ (and then $\mu(a)\neq x$, because $\mu(x)\neq c$);
2) or $\mu(y)=a$ and, correspondingly, $\mu(a)=x$ (and $y\neq c$).
Let us first consider case 2).

Since $\mu(b)\neq x$, in this case, the triple $(c,b,x)$ is not blocking only
if $\mu(c)=b$ (see red edges in Fig. 1). Here we again get two alternatives:
a) $y\neq d$. Then $R_{\mu}(d)>2$. Since $R_{\mu}(b)>1$, this means that the
triple $(d,b,x)$ is blocking.
b) $y=d$. Then $\mu(y)=a\in\\{a,b,t\\}$, which was to be proved.

$c$$d$$e$$a$$t$$x$$b$$s$$r^{\prime}_{x}$$r_{x}^{\prime}+1$33 Fig. 2: The part
of the graph defining preferences that is considered in Lemma 1. Red lines
represent edges that correspond to pairs that belong to families in the
matching $\mu$ in case I.1.b

Thus, in case 2), the assertion of the lemma is valid. Let us consider case
1), i.e., assume that $\mu(c)=a$ and $\mu(a)\neq x$. In this case, two
alternatives are possible:
a) $y\neq d$. Since $R_{\mu}(a)>1$ and $R_{\mu}(d)>1$, the triple $(d,a,x)$ is
blocking.
b) $y=d$. In this case, the triple $(d,a,e)$ is not blocking only if
$\mu(a)=e$. Then the matching $\mu$ contains the family $(c,a,e)$, while
$\mu^{-1}(x)\neq b$ (otherwise this matching would have contained $(d,b,x)$).
Then $R_{\mu}(b)>2$ and $R_{\mu}(d)>2$, consequently, the triple $(d,b,s)$ is
not blocking only if $\mu(b)=s$ (here $\mu(s)\neq d$, because $d=y$), see Fig.
2. Since $\mu^{-1}(x)\neq t$ (otherwise the matching would have contained the
family $(d,t,x)$), in the case under consideration, $R_{\mu}(d)>3$. Moreover,
$\mu(t)\not\in\\{e,s\\}$, therefore, $R_{\mu}(t)>2$. We can analogously prove
that $R_{\mu}(s)>1$. But then the triple $(d,t,s)$ is blocking.

II) It remain to study the case, when the initial assumption that $\mu(x)\neq
c$ is false.

Let $\mu(x)=c$. Assume that $\mu(c)\not\in\\{a,b\\}$. Then $R_{\mu}(c)>2$,
$R_{\mu}(a)>1$, and $R_{\mu}(b)>1$. Consider the triple $(c,a,e)$. It is not
blocking only if $\mu(a)=e$. But then blocking is the triple $(c,b,e)$.
$\square$

Note that in the statement of Lemma 1, we can use a subgraph with the same
number of vertices, but with a lesser number of edges (for example, the graph
without the edge $(t,e)$, but with the edge $(t,s)$ of rank 1). We need Lemma
1 (in the proved form) for verifying Lemma 2.

Let the symbol $H$ stand for the graph of the 3DSMI problem. For any vertex
$v\in V(H)$ we put $\rho_{H}(v)=\max\limits_{(v,v^{\prime})\in
E(H)}r(v,v^{\prime}).$.

###### Theorem 1

Let $H$ be the graph of the 3DSMI problem of size $n$. Let us construct from
it the graph $G$ of the 3DSM problem as follows. The graph $H$ is a subgraph
of the graph $G$ (with the same ranks of edges). To each vertex $x$ of the
graph $H$ we “attach” the corresponding copy of the graph $H_{x}$ shown in
Fig. 1 with vertices $a_{x},b_{x},c_{x},d_{x},e_{x},s_{x},t_{x}\not\in V(H)$
(all subgraphs $H_{x}$ are pairwise disjoint). Moreover, let the value
$r^{\prime}_{x}$ in formulas (1) equal $\rho_{H}(x)+1$. Let us define ranks of
the rest edges of the graph $G$ of the 3DSM problem arbitrarily. Then, if the
graph $H$ of the 3DSMI problem has no stable matching, then neither does the
constructed graph $G$ of the 3DSM problem. Here the size of the considered
3DSM problem equals $8n$.

Proof (cf. proofs of lemmas 2 and 6 in [10]): Assume the contrary, i.e., for
the 3DSM problem with the graph $G$ there exists a stable matching $\mu_{G}$.
Then for any vertex $x\in V(H)$ the assertion of Lemma 1 is valid. In other
words, either alternative A takes place, while $y=\mu_{G}(x)\not\in V(H)$ and
$\mu_{G}(y)\not\in V(H)$, or alternative B takes place, while $(x,y)\in E(H)$
and $(y,\mu_{G}(y))\in E(H)$ (the violation of the latter correlation, i.e.,
the fulfillment of alternative A for the vertex $y$, leads to the
contradiction: $x=\mu_{G}(\mu_{G}(y))\in\\{a_{y},b_{y},t_{y}\\}\not\in V(H)$).

Note that if alternative A takes place, then

$\rho_{H}(x)<R_{\mu_{G}}(x).$ (2)

Let us associate the matching $\mu_{G}$ with the matching $\mu_{H}$ of the
3DSMI problem with the graph $H$. Assume that in the case of alternative A,
$\mu_{H}(x)=x$ (i.e., the agent $x$ remains single). In the case of
alternative B, we put $\mu_{H}(x)=\mu_{G}(x)$,
$\mu_{H}(\mu_{H}(x))=\mu_{G}(\mu_{G}(x))$.

By condition, the matching $\mu_{H}$ is not stable. Therefore, there exists a
blocking triple $(v,v^{\prime},v^{\prime\prime})$, where
$v,v^{\prime},v^{\prime\prime}\in V(H)$. The blocking property means that each
vertex in this triple satisfies the inequality for ranks, for example, for the
vertex $v$ this inequality takes the form $r(v,v^{\prime})<R_{\mu_{H}}(v)$. If
$\mu_{H}(v)\neq v$, then $\mu_{G}(v)=\mu_{H}(v)$ and, consequently,
$r(v,v^{\prime})<R_{\mu_{G}}(v)$. But if $\mu_{H}(v)=v$, then
$r(v,v^{\prime})<R_{\mu_{G}}(v)$ due to inequality (2). In any case, the same
triple $(v,v^{\prime},v^{\prime\prime})$ is blocking for the matching
$\mu_{G}$ in the 3DSM problem, i.e., $r(v,v^{\prime})<R_{\mu_{G}}(v)$,
analogously, $r(v^{\prime},v^{\prime\prime})<R_{\mu_{G}}(v^{\prime})$ and
$r(v^{\prime\prime},v)<R_{\mu_{G}}(v^{\prime\prime})$.

Let us calculate the size of the 3DSM problem. Each triple of vertices of the
graph $H$, which are associated with different genders, corresponds to three
various subgraphs in the form shown in Fig. 1. In one of the mentioned cases,
the vertex $x$ and all vertices of the same color (as $x$) shown in Fig. 1 are
colored green (therefore, we get three green, three blue, and two red
vertices). In the other case, the vertex $x$ is colored blue, while green are
only vertices $c_{x}$ and $d_{x}$ (therefore, we get two green, three blue,
and three red vertices) and, finally, in the remaining case, the subgraph
contains three green, two blue, and three red vertices. Therefore, the number
of vertices that correspond to each gender increases by 8 times. $\square$

##  4 An example of an unsolvable 3DSM problem of size 18

As was mentioned above, Theorem 1, along with the result obtained in the paper
[11] (see Lemma 3 in the Appendix), allows one to construct instances of the
3DSM problem of size 24 with no stable matching. In this section, we reduce
the size of such a counterexample to 18. As a base we use a certain
modification of Lemma 1.

$c$$d$$e$$f$$a$$t$$x$$b$$s$$z$33433 Fig. 3: The part of the graph that is
considered in Lemma 2. The rank of bold edges equals 1. Dashed lines represent
edges of rank 2. Edges $(x,c)$ and $(x,d)$ have ranks $r^{\prime}_{x}$ and
$r^{\prime}_{x}+1$, edges $(z,c)$ and $(z,d)$ do ranks $r^{\prime}_{z}$ and
$r^{\prime}_{z}+1$, edges $(c,f)$, $(b,e),(a,e)$, and $(d,t)$ have rank $3$,
$r(b,s)=4$

Recall that in Lemma 1 one considers a subgraph, whose copy we “attach” to
each vertex of the graph $H$ that defines the 3DSMI problem with no stable
matching. In Lemma 2, one considers a subgraph, which in certain cases can
have “two attachments” to two vertices of such a graph $H$. The “cost” of this
effect is the supplement of the subgraph with the vertex $f$ (apart from the
“attached” vertex $z$); cf. Fig. 1 and Fig. 3.

###### Lemma 2

Assume that in the 3DSM problem the subgraph of the graph $G$ takes the
following form (see Fig. 3) with indicated ranks of edges. Let $\mu$ be a
stable matching. Then
A) either $\mu(x)\in\\{c,d\\}$ and $\mu(\mu(x))\in\\{a,b,f,t\\}$;
B) or $R_{\mu}(x)<r^{\prime}_{x}$ and then a) either $\mu(z)\in\\{c,d\\}$ and
$\mu(\mu(z))\in\\{a,b,f,t\\}$,
ttt b) or $R_{\mu}(z)<r^{\prime}_{z}$.

Proof: Let us prove that

$\text{either alternative~{}A) takes place or }R_{\mu}(x)<r^{\prime}_{x}.$ (3)

Let us first prove that the inequality $\mu(x)\neq c$ guarantees the
fulfillment of these conditions.

Assume the contrary. Repeating verbatim part I of the proof of Lemma 1, we
note that if $\mu(c)=a$ (i.e., $\mu(x)\neq c$), then it suffices to consider
the case, when $\mu(x)=d$ (i.e., case I.2.b in Lemma 1). In this case, the
triple $(d,a,z)$ is blocking, only if the matching $\mu$ does not contain the
family $(c,a,z)$. Moreover, by our assumption, the triple $(d,b,x)$ cannot
belong to this matching. Then $R_{\mu}(b)>2$ and $R_{\mu}(d)>2$. Consequently,
the triple $(d,b,e)$ is blocking, only if the equality $\mu(b)=e$ is false.
Since $\mu^{-1}(x)\neq t$ (otherwise the matching would have contained the
family $(d,t,x)$, we recall that $\mu(x)=d$), in the case under consideration,
$R_{\mu}(d)>3$. In addition, $\mu(e)\not\in\\{c,d\\}$, therefore,
$R_{\mu}(e)>2$. Analogously, $R_{\mu}(t)>1$. But then the triple $(d,t,e)$ is
blocking (see Fig. 4).

$c$$d$$e$$f$$a$$t$$x$$b$$s$$z$33433 Fig. 4: The part of the graph defining
preferences that is considered in Lemma 2. Red lines represent the edges
corresponding to pairs which belong to families in the matching $\mu$ in case
I.2.b in Lemma 1.

For verifying proposition (3), it remains to prove that in the case, when
$\mu(x)=c$, we get the inclusion $\mu(c)\in\\{a,b,f\\}$. Assume that
$\mu(x)=c$ and $\mu(c)\not\in\\{a,b,f\\}$. Then $R_{\mu}(c)>3$ and
$R_{\mu}(e)>1$. Both triples $(a,e,c)$ and $(b,e,c)$ are not blocking, only if
$\\{\mu^{-1}(e),\mu^{-1}(z)\\}=\\{a,b\\}$. But then blocking is the triple
$(f,e,c)$.

We have proved the validity of proposition (3). Let us now prove that in the
case, when $R_{\mu}(x)<r^{\prime}_{x}$,

one of alternatives a) or b) mentioned in assumptions of Lemma 2 takes place.
(4)

If $R_{\mu}(x)<r^{\prime}_{x}$, then it makes no sense to consider families
that contain the vertex $x$ and other vertices shown in Fig. 3. Therefore, we
can delete from the graph the vertex $x$, together with the edges of the graph
that are incident to it. At the same time, we can change ranks of all edges
that go out of the vertices that previously represented initial points of
edges with the endpoint $x$, preserving the order of ranks. In other words, we
can reduce ranks of all edges, which go out of vertices $a$ and $b$, by one.
This operation does not affect the validity of proposition (4).

If we use the symbol $x$ instead of $z$, then the graph in the modified in
such away Fig. 3 will contain the subgraph shown in Fig. 1. According to Lemma
1, either alternative a) or b) takes place. $\square$

We need a concrete instance of the 3DSMI problem of size 3 with no stable
matching proposed in [11, Theorem 2]. The graph $H$ of the considered problem
is shown in Fig. 5. The 3DSMI problem that corresponds to this graph has no
stable matching. See Appendix for the proof of this fact.

501234678 Fig. 5: The graph of the 3DSMI problem of size 3 with no stable
matching. For convenience, we denote vertices of the graph $H$ by numbers $v$,
$v=0,1,\ldots,8$. The value $v\bmod 3$ defines the gender that corresponds to
the vertex $v$. The rank of bold edges equals 1. Dashed lines represent the
edges of rank 2. The rank of dotted edges equals 3.

###### Theorem 2

Let the graph $G$ that corresponds to the 3DSM problem contain the subgraph
$H$ shown in Fig. 5. Assume that in all considered below subgraphs of the
graph $G$, which are identical to graphs mentioned in lemmas 1 and 2, the
value $r^{\prime}_{s}$ in the corresponding copies coincides with
$\rho_{H}(s)+1$, where $s$ is some vertex (we specify its number below, in
assumptions of lemmas 1 and 2 we mean either the vertex $x$ or $z$). Thus, the
graph $G$ contains three nonintersecting subgraphs $H_{0}$, $H_{2}$, and
$H_{7}$, which are identical to the graph mentioned in lemma 1. The vertex $x$
in these graphs represents, correspondingly, vertices 0, 2, and 7 of the graph
$H$, and subgraphs $H_{0}$, $H_{2}$, and $H_{7}$ have no other common vertices
with the graph $H$. Moreover, the graph $G$ contains three nonintersecting
subgraphs $H_{1,4}$, $H_{3,6}$, and $H_{5,8}$, which are identical to
subgraphs mentioned in Lemma 2. The vertex $x$ in these subgraphs represents,
correspondingly, vertices 1, 3 and 5, while the vertex $z$ does vertices 4, 6,
and 8; subgraphs $H_{1,4}$, $H_{3,6}$, and $H_{5,8}$ have no other common
vertices with the graph $H$ and the graphs $H_{0}$, $H_{2}$, $H_{7}$. Assume
that the graph $G$ has no vertices different from those considered above,
namely,

$V(G)=\bigcup_{v\in\\{0,2,7\\}}V(H_{v})\bigcup_{(v,v^{\prime})\in\\{(1,4),(3,6),(5,8)\\}}V(H_{v,v^{\prime}}).$

We can arbitrarily define ranks of edges of the graph $G$ that were not
mentioned above. Then the 3DSM problem defined by the graph $G$ has no stable
matching.

Proof: Assume the contrary, namely, assume that for the graph $G$ there exists
a stable matching $\mu_{G}$. Assume also that for vertex 1 alternative A
mentioned in Lemma 2 takes place, i.e., $\mu_{G}(1)\in\\{c_{1},d_{1}\\}$ and
$\mu_{G}^{-1}(1)\in\\{a_{1},b_{1},f_{1},t_{1}\\}$. The triple $(0,1,5)$ is not
blocking only if $\mu_{G}(5)=0$. Then for vertex 0, alternative A in Lemma 1
is impossible; therefore, $\mu_{G}(0)=7$, i.e., the matching $\mu_{G}$
contains the family $(7,5,0)$. Therefore, for vertex 7 neither alternative A
nor B in Lemma 1 is possible. The obtained contradiction means that for vertex
1 alternative B in Lemma 2 takes place.

Assume now that vertex 3 satisfies alternative A in Lemma 2. Then the triple
$(1,2,3)$ is not blocking, only if $\mu_{G}(1)=2$. But then vertex 2 can
satisfy neither alternative A no B in Lemma 1. The obtained contradiction
means that vertex 3 satisfies alternative B in Lemma 2.

Assume that vertex 5 satisfies alternative A in Lemma 2. Then the triple
$(3,4,5)$ is not blocking, only if $\mu_{G}(3)=4$. But the equality
$\mu_{G}(1)=5$ is impossible for vertex 1, because
$\mu_{G}^{-1}(5)\in\\{a_{5},b_{5},f_{5},t_{5}\\}$. Therefore, $\mu_{G}(1)=2$,
and vertex 2 satisfies alternative B in Lemma 1. Then the matching $\mu_{G}$
contains the family $(1,2,3)$, which contradicts the assumption that
$\mu_{G}(3)=4$. Consequently, vertex 5 also satisfies alternative B in Lemma
2.

Therefore, for all vertices $x\in V(H)$, one of two alternatives takes place:
either $\mu_{G}(x)\not\in V(H)$ and $\mu_{G}^{-1}(x)\not\in V(H)$; or
$(x,\mu_{G}(x))\in E(H)$ (in particular, the latter case takes place for
vertices $x=1,3,5$). Evidently, in the case of the second alternative, for
$y=\mu_{G}(x)$ we also get the inclusion $(y,\mu_{G}(y))\in E(H)$ (cf. the
same correlation in the proof of Theorem 1).

Let us associate the stable matching $\mu_{G}$ with the matching $\mu_{H}$ in
the 3DSMI problem with the graph $H$ in the same way as in the proof of
Theorem 1. Namely, in the case, when the first alternative is fulfilled for
$x\in V(H)$, we put $\mu_{H}(x)=x$; in the case of the second one, we do
$\mu_{H}(x)=\mu_{G}(x)$, $\mu_{H}(\mu_{H}(x))=\mu_{G}(\mu_{G}(x))$.

Since for the graph $H$ no stable matching exists, for the matching $\mu_{H}$,
one can find a blocking triple $(v,v^{\prime},v^{\prime\prime})$, where
$v,v^{\prime},v^{\prime\prime}\in V(H)$. Similarly to the proof of Theorem 1,
we conclude that the same triple $(v,v^{\prime},v^{\prime\prime})$ is blocking
in the matching $\mu_{G}$ in the 3DSM problem. $\square$

Let us calculate the size of the 3DSM problem mentioned in Theorem 2. Each of
subgraphs $H_{0}$, $H_{2}$, and $H_{7}$, has, in total, 8 vertices that
correspond to various genders, while each subgraph among $H_{1,4}$, $H_{3,6}$,
and $H_{5,8}$ does 10 ones. Therefore, the graph $G$ mentioned in Theorem 2
defines the 3DSM problem of size 18.

##  5 Concluding remarks and open problems

In this paper, we study the correspondence between 3DSMI and 3DSM problems
with no stable matching. This idea belongs to C.K. Lam and C.G. Plaxton; they
were first to construct an instance of the 3DSM problem of size $n=90$ with no
stable matching [10]. We use the graph theory language, which is more
demonstrative. Finally, we have succeeded in reducing the size of the
counterexample for the 3DSM problem by 5 times in comparison with the initial
value of $n=18$.

Earlier K. Pashkovich and L. Poirrier ([13]) have proved that in any 3DSM
problem of size $n\leqslant 5$ there exists a stable matching.

Therefore, the minimum size of the 3DSM problem with no stable matching is not
less than 6 and not greater than 18. Its exact value is not known yet.
Possibly, computer experiments, including minimal unsolvable 3DSM problems
stated by us earlier in [11], can be helpful in answering this question.

Note that in experiments with randomly generated 3DSM problems, whose size $n$
is close to 20 (with the independent uniform distribution of all possible
permutations in preference lists), heuristic algorithms always succeed in
finding a stable matching within a relatively short time. Examples of 3DSM
problems with no stable matching seem to be rather rare. The following
proposition proved by B. Pittel indirectly indicates this fact: the mean value
of stable matchings grows faster than $n^{2}\ln^{2}(n)$ as $n\to\infty$
([14]). So the following question seems to be interesting: does the fraction
of unsolvable 3DSM problems tend to zero as $n\to\infty$?

In the classical 2DSM problem, stable matchings always exist; moreover,
according to the paper [9] by J.H. Conway and D.E. Knuth, the introduction of
the natural order relation makes them form a distribution lattice. Moreover,
as was proved later [2, 7], the correspondence between distribution lattices
and stable matchings is biunique. The authors of the mentioned papers have
also explicitly described a partially ordered set (POSET), whose ideals define
a lattice of stable matchings. Elements of this POSET represent the so-called
rotations (i.e., cycles that include alternately men and women) defined by
preference matrices of the 2DSM problem [8].

We are interested in a “correct” statement of this problem in the case of
three genders (3G), where stable matchings would have the same elegant
combinatorial structure. As is well known, stable matchings always exist in
the 3GSM problem with V.I. Danilov lexicographically acyclic preferences [3].
Let us give another statement of the problem, in which stable matchings always
exist. To this end, we define 6 matrices of pairwise preferences between
genders; we say that a matching is stable, if there exists no triple, which
does not belong to this matching, where pairwise relations of each pair in a
triple “do not get worse”. However, the combinatorial structure of stable
matchings in the “correct” three-gender problem needs to be studied yet.

## References

  * [1] P. Biró, E. McDermid, Three-sided stable matchings with cyclic preferences, Algorithmica 58 (2010) 5–18.
  * [2] C. Blair, Every finite distributive lattice is a set of stable matchings, Journal of Combinatorial Theory, Series A, 37 (3) (1984) 353–356
  * [3] V.I. Danilov, Existence of stable matchings in some three-sided systems. Mathematical Social Sciences. 46(2) (2003) 145–148.
  * [4] E. Boros, V. Gurvich, S. Jaslar, D. Krasner, Stable matchings in three-sided systems with cyclic preferences. Discrete Math. 289(1’3) (2004) 1–10.
  * [5] K. Eriksson, J. Söstrand, P. Strimling, Three-dimensional stable matching with cyclic preferences, Math. Soc. Sci. 52 (2006) 77–87.
  * [6] D. Gale, L.S. Shapley, College admissions and the stability of marriage, Am. Math. Mon. 69 (1962), 9–15.
  * [7] D. Gusfield, R.W. Irving, P. Leather and M. Saks, Every finite distributive lattice is a set of stable matchings for a small stable marriage instance. J. Combinatorial Theory A44 (1987) 304-309.
  * [8] R.W. Irving, P. Leather, The complexity of counting stable marriages. SIAM Journal on Computing, 15(3) (1986) 655–667
  * [9] D.E. Knuth, Stable marriage and its relation to other combinatorial problems: an introduction to the mathematical analysis of algorithms, in: CRM Proceedings and Lecture Notes, 1996.
  * [10] C.-K. Lam, C.G. Plaxton, On the Existence of Three-Dimensional Stable Matchings with Cyclic Preferences, in: Lecture Notes in Computer Science, 11801, 2019. Algorithmic Game Theory, 329–342.
  * [11] E.Yu. Lerner, R.E. Lerner, Minimal instances with no weakly stable matching for three-sided problem with cyclic incomplete preferences, arXiv.org. 2021. arXiv:2101.08223 [math.CO].
  * [12] D.F. Manlove, Algorithmics of matching under preferences, Theor. Comput. Sci. World Scientific, 2013.
  * [13] K. Pashkovich, L. Poirrier, Three-dimensional stable matching with cyclic preferences, Optimization Letters (2020).
  * [14] B. Pittel, On random stable matchings: Cyclic ones with strict preferences and two-sided ones with partially ordered preferences, Advances in Appl. Math. 120 (2020) 1–27.

## Appendix. An instance of the 3DSMI problem of size $n=3$ with no stable
matching

######  Lemma 3 ([11], Theorem 2)

The 3DSMI problem with the graph $G$ shown in Fig. 5 has no stable matching.

Proof: There exist 7 families that form matchings in this problem, namely,
$(0,1,5)$, $(0,7,8)$, $(1,2,3)$, $(1,5,3)$, $(2,3,4)$, $(3,4,5)$, and
$(4,8,6)$.

Recall that a matching $\mu$ in the 3DSMI problem defined by the graph $G$ is
said to be complementable, if there exists a triple of vertices
$(v,v^{\prime},v^{\prime\prime})$ such that $\mu(v)=v$,
$\mu(v^{\prime})=v^{\prime}$, $\mu(v^{\prime\prime})=v^{\prime\prime}$, and
$\\{(v,v^{\prime}),(v^{\prime},v^{\prime\prime})(v^{\prime\prime},v)\\}\subseteq
E(G)$.

Evidently, any complementable matching is not stable. For the graph shown in
Fig. 5 there exists 8 noncomplementable matchings. Below we give their
complete list together with blocking triples:
1) $\\{(0,1,5),(2,3,4)\\}$, the blocking triple is $(4,8,6)$;
2) $\\{(0,1,5),(4,8,6)\\}$, the blocking triple is $(1,2,3)$;
3) $\\{(0,7,8),(1,2,3)\\}$, the blocking triple is $(3,4,5)$;
4) $\\{(0,7,8),(1,5,3)\\}$, the blocking triple is $(2,3,4)$;
5) $\\{(0,7,8),(2,3,4)\\}$, the blocking triple is $(0,1,5)$;
6) $\\{(0,7,8),(3,4,5)\\}$, the blocking triple is $(0,1,5)$;
7) $\\{(1,2,3),(4,8,6)\\}$, the blocking triple is $(0,7,8)$ or $(3,4,5)$;
8) $\\{(1,5,3),(4,8,6)\\}$, the blocking triple is $(0,7,8)$. $\square$