# Finite Elasticity of the Vertex Model and its Role in Rigidity of Curved
Cellular Tissues.

Arthur Hernandez<EMAIL_ADDRESS>Department of Physics, University
of California Santa Barbara, Santa Barbara, CA 93106  Michael F. Staddon
Center for Systems Biology Dresden, Dresden, Germany  Max Planck Institute for
the Physics of Complex Systems, Dresden, Germany  Max Planck Institute of
Molecular Cell Biology and Genetics, Dresden, Germany  Michael Moshe  Racah
Institute of Physics, The Hebrew University of Jerusalem, Jerusalem, Israel
91904  M. Cristina Marchetti  Department of Physics, University of California
Santa Barbara, Santa Barbara, CA 93106

###### Abstract

Using a mean field approach and simulation, we study the non-linear mechanical
response of the vertex model (VM) of biological tissue under compression and
dilation. The VM is known to exhibit a transition between rigid and fluid-
like, or floppy, states driven by geometric incompatibility. Target perimeter
and area set a target shape which may not be geometrically achievable, thereby
engendering frustration. Previously, an asymmetry in the linear elastic
response was identified at the rigidity transition between compression and
dilation. Here we show and characterize how the asymmetry extends away from
the transition point for finite strains. Under finite compression, an
initially solid VM can totally relax perimeter tension, and thereby have
reduced bulk and shear modulus. Conversely, an initially floppy VM under
dilation can rigidify and have a higher bulk and shear modulus. These
observations imply that re-scaling of cell area shifts the transition between
rigid and floppy states. Based on this insight, we calculate the re-scaling of
cell area engendered by intrinsic curvature and write a prediction for the
rigidity transition in the presence of curvature. The shift of the rigidity
transition in the presence of curvature for the VM provides a new metric for
predicting tissue rigidity from image data for curved tissues in a manner
analogous to the flat case.

††preprint: APS/123-QED

##  I Introduction

Understanding the emergence of form in organ development presents a major
challenge to current continuum physics modeling of biological tissues. Unlike
passive materials, certain tissues may tune their mechanical response to
applied strains and forces by modifying cell shape and thereby be rigid or
floppy. In particular, cell shape as characterized by the shape index
$s\equiv\frac{P}{\sqrt{A}}$ has been found to serve as a signal for the onset
of a solid-liquid transition at constant density in epithelial tissues [1].

One widely studied model of epithelial tissues is the vertex model (VM) which
models the epithelium as a collection of vertices and edges in the 2D plane,
and reduces the tissue’s structural data to a polygonal tiling with possible
edge tension. Unlike conventional spring network models which penalize
deviations away from each edge length, the VM instead only sets a target cell
area due to 3D bulk tissue incompressibility, along with terms capturing cell-
cell edge adhesion and cell contractility[2], which constrain the cell’s
perimeter. Thus the VM is less constrained than a spring network, e.g.
crystalline solids, and naturally engenders zeros modes for any polygonal
tiling.[3, 4, 2].

What’s more these zero modes exist at the level of a single cell. For example,
assuming all cell edges have identical adhesion and contractility, the VM
energy reduces to penalizing harmonic deviations away from a target area
$A_{0}$ and target perimeter $P_{0}$. Thus in the VM each cell has 2 shape
constraints, but a general polygon has at least 3 degrees of freedom such as
is the case for triangles [3].

Based on constraint counting it seems the VM can never support a rigid state.
Nonetheless the VM exhibits a rigidity transition between solid and floppy
states tuned by the target shape index $s_{0}=\frac{P_{0}}{\sqrt{A_{0}}}$
about a critical point $s_{0}^{*}$ [5]. The transition is due to a geometric
constraint set by the isoperimetric inequality which gives a lower bound on
the ratio of $\frac{P}{\sqrt{A}}$ for n-gons admissible on the plane [6].

$\displaystyle\frac{P}{\sqrt{A}}\geq s_{0}^{*}(n)$ (1)

Where $s_{0}^{*}=\sqrt{4n\tan\left(\frac{\pi}{n}\right)}$ is the isoperimetric
quotient. The lower bound sets an incompatible regime
$\frac{P_{0}}{\sqrt{A_{0}}}<s_{0}^{{*}}$ where polygons cannot simultaneously
achieve $A_{0}$ and $P_{0}$, and a compatible regime
$\frac{P_{0}}{\sqrt{A_{0}}}\geq s_{0}^{{*}}$ where polygons may achieve both
$A_{0}$ and $P_{0}$. This geometric constraint on shape indicates that
rigidity stems from self-tension due to geometrically incompatibility.

For VM simulations consisting of ordered tilings (triangles, squares,
hexagons) in the plane the rigidity transition occurs exactly at
$s_{0}^{*}(n)$ [3]. Whereas for disordered VM simulations the transition
occurs at approximately $s_{0}^{*}(5)$ [5].

The linear response of the VM to mechanical deformations is well studied [7,
3, 8, 9] but the non-linear response relatively less so. Recent work showed
that the VM exhibits shear-thickening in the compatible regime [4]. In the
same vein, this paper presents a careful study of the non-linear elasticity of
the VM under finite dilation and compression via a mean-field approach and
simulation. In previous work by the authors, the onset of compatibility in the
VM at $s_{0}^{*}$ showed anomalous elasticity as reflected by an asymmetric
bulk modulus under dilation and compression, as well as coupling between
stretching and shear modes [3].

In this article, we show that the asymmetry of the bulk modulus extends away
from $s_{0}^{*}$ under finite compression and dilation. In particular, the VM
exhibits a dilation-hardening for compatible tissues and a compression-
softening in incompatible tissues for finite critical strain. The hardening
(softening) nonlinear response to dilation (compression) is reflected by a
jump (drop) discontinuity of the bulk modulus and is associated with the
sudden lifting (onset) of zero-modes.

The mechanism of a tissue fine-tuning their rigidity in response to areal re-
scaling relates to several biological processes such as tissue growth,
shrinkage, applied deformations, and in particular spontaneous generation of
Gaussian (intrinsic) curvature. Based on insight from the planar 2D non-linear
elasticity, we use of mean field theory to predict how local
compression/dilation due to intrinsic curvature shifts the transition point
between rigid and floppy states.

The organization of the paper is as follows: in section II we define our mean-
field approach of the VM which models 2D tissue elasticity at the single cell
level. Section III outlines the calculation of the non-linear bulk modulus of
the mean-field model and discusses simulation results. In section IV we
present a Landau energy argument to elucidate the connection between the
asymmetry of the bulk modulus and how a finite critical strain controls the
onset/lifting of zero-modes. Section V uses our results from the mean field
theory to predict the effective critical shape index for cells on a curved
surface. We compare our prediction for the rigidity transition in the presence
of curvature to simulations done by [10]. We conclude with discussion in
section VI.

##  II Mean-Field Theory of ordered vertex model

Our mean field theory for the VM is defined as a uniform regular 2D tiling
with all cells responding identically to applied deformations. The tissue
energy is the sum of all individual cells, and therefore our mean field
approach reduces the VM to a single polygonal cell. All bulk tissue
properties, such as elastic moduli, are defined at the single cell level
response. Details of the mean field model are given in appendix C, and a
thorough study by the authors is in [3]. For concreteness, our simulations and
mean field theory are for hexagonal cells unless stated otherwise. All results
hold analogously for other polygons.

The tissue energy per cell is rescaled by $\kappa_{A}A_{0}^{2}$ such that the
energy only has two dimensionless free parameters. The VM energy per cell is
cast as

$\displaystyle
E=\frac{1}{2}\left(a-1\right)^{2}+\frac{r}{2}\left(p-s_{0}\right)^{2}$ (2)

where $r\equiv\frac{\kappa_{P}}{\kappa_{A}A_{0}}$ is the rigidity ratio, and
$s_{0}\equiv\frac{P_{0}}{\sqrt{A_{0}}}$ the target shape index, and $a$ and
$p$ are the scaled area and perimeter.

To parameterize cell shape degrees of freedom, we work with a Cartesian
coordinate system $(X,Y)$ encompassing the cell with $Y$ along the height, and
$X$ along the width. The area and perimeter of a cell are purely geometric
objects, and their shape change under various deformations can be computed
given a transformation law.

The externally imposed dilation and compression are implemented via an overall
re-scaling of the cell’s height $h$ and width $w$ via the transformation
$w\rightarrow w(1+\epsilon)$, $h\rightarrow h(1+\epsilon)$, where
$\epsilon\in(-1,1)$. In response to the strain, the cell may also
spontaneously shear while maintaining the imposed rescaled area, as shown in
Fig 1. This ”tilt” response is parametrized via a self-shear deformation, with
$w\rightarrow w+th$, $h\rightarrow h$, where
$t(\theta)\equiv\tan\left({\theta}\right)$, and represents the shape
degeneracy of cells.

Figure 1: Under compression/dilation, the cell may respond via an affine self-
shear transformation where it tilts either right or left. Only parallel edges
change length in the same manner, whereas neighboring edges change length
differently.

In addition, hexagonal cells can also respond via non-affine transformations,
and what’s more these non-affine responses yield softer shear and Young’s
moduli in the incompatible regime [11]. In this study we preclude non-affine
pathways for cell response as our previous work showed this approximation well
captures the response under imposed isotropic compression/dilation [3].

The deformed energy of an isotropically dilated or compressed cell is then
given by

$\displaystyle E(\epsilon,\theta;s_{0},r)=$
$\displaystyle\frac{1}{2}\left[\ell^{2}(1+\epsilon)^{2}-1\right]^{2}+\frac{r}{2}\left[p(\epsilon,t(\theta),\ell)-s_{0}\right]^{2}\;,$
(3)

where $\ell(s_{0},r)$ is the rescaled ground state characteristic cell size
(see appendix C for details) and the deformed perimeter is

$\displaystyle p(\epsilon,\theta,\ell)$
$\displaystyle=\frac{\sqrt{2}\ell}{3^{3/4}}(1+\epsilon)\bigg{(}2\left(1+t\right)+\sqrt{1+(t(\theta)-\sqrt{3})^{2}}$
$\displaystyle+\sqrt{1+(t(\theta)+\sqrt{3})^{2}}\bigg{)}\;.$ (4)

If we set $\epsilon=0$ and minimize with respect to $\theta$ we recover the
previous ground states results in [3] where the ground state energy is gapped
for $s_{0}<s_{0}^{*}$ and vanishes for $s_{0}\geq s_{0}^{*}$, populated by a
manifold of degenerate shapes parametrized by $\theta$. To study the response
to $\epsilon\neq 0$ we minimize $\theta$ as a function of applied strain in a
manner analogous to our study of the linear response [3]. Formally, the
energetic response is given by

$\displaystyle E=\min_{\theta}E(\theta;\epsilon,s_{0},r)$ (5)

Because height and width are fixed by dilation/compression the energy
minimization is 1D and corresponds to solving,

$\displaystyle\frac{\partial
E}{\partial\theta}\bigg{|}_{s,r,\epsilon}=\left(p-s_{0}\right)\frac{\partial
p}{\partial\theta}\bigg{|}_{s_{0},r,\epsilon}=0$ (6)

Which has two solutions: either a cell utilizes shape degeneracy via $\theta$
so that the perimeter accommodates both dilation/compression and target shape
index $s_{0}$, or the perimeter is totally set by dilation/compression with no
tilt response. The relevant energy minimizing solution is a function of
$s_{0},r$ and $\epsilon$.

##  III Nonlinear elasticity in the Euclidean Plane

The non-linear response under finite dilation and compression is characterized
by the bulk modulus, defined as

$\displaystyle
K=\frac{1}{2a_{\text{cell}}}\left(\frac{\partial^{2}}{\partial\epsilon^{2}}\min_{\theta_{\text{min}}}E(\epsilon,\theta;s_{0},r)\right)_{s_{0},r,\epsilon}$
(7)

where $a_{\text{cell}}=\frac{3\sqrt{3}}{2}\ell^{2}$ is the rescaled cell area.
Evaluating Eq. 7 at $\epsilon=0$ yields the linear response, whereas a finite
$\epsilon$ gives the non-linear response under finite strains.

It is crucial that the minimization occurs before differentiation because the
self-shear response $\theta$ is implicitly dependent on $\epsilon$ via Eq. 6.

In the incompatible rigid state, ($s_{0}<s_{0}^{*}$), we find that the mean
field model and simulation exhibit a discontinuous drop in the bulk modulus at
a critical compression. The discontinuity occurs due to a spontaneous self-
shear of the cell which allows perimeter tension to vanish. Conversely, under
dilation the bulks modulus remains continuous as a function of strain, see
Fig. 2 Plot A. Larger choices of $r$ will shift the critical strain up,
reflecting how a higher perimeter tension may support higher compression
before giving way to spontaneous self-shear.

In the compatible floppy state, ($s_{0}>s_{0}^{*}$), the bulk modulus is
continuous under any finite compression but discontinuous jump at some
critical dilation, see Fig. 2 plot B. At sufficient dilation, the zero-modes
of the degenerate ground state are ”exhausted”, resulting in a frustrated and
thereby rigid state. Unlike the incompatible state, the critical dilation is
insensitive to $r$.

At the transition, $s_{0}^{*}$, both dilation hardening and compression
softening are present for arbitrarily small strains, and reflect an asymmetry
of the response to area re-scaling. For various dilation/compression strain
magnitudes we plot the difference between dilation and compression bulk
modulus $\Delta K$ vs. $s_{0}$ as a measure of asymmetry of the response in
Fig. 2 plot C. Around the critical shape index, $s_{0}^{*}$, the asymmetry
extends continuously away from the critical point for even modest strain
values ($>0.002$). Note the curve $\Delta K$ vs. $s_{0}$ is also not left-
right symmetric along the $s_{0}=s_{0}*$ axis; this is due to the critical
strain only depending on the choice of $r$ in the incompatible state but not
in the compatible state.

The origin of the bulk modulus discontinuity can be partial elucidated by
writing out Eq.8 explicitly

$\displaystyle
K=\frac{1}{2a_{\text{cell}}}\left(a^{2}+r\big{(}p(\theta_{min})-s_{0}\big{)}\frac{\partial^{2}p(\theta_{min})}{\partial\epsilon^{2}}+r\left(\frac{\partial
p(\theta_{min})}{\partial\epsilon}\right)^{2}\right).$ (8)

In the incompatible state, if cells can accommodate target perimeter and
imposed strains simultaneously, the perimeter tension,
$r\big{(}p(\theta_{min})-s_{0}\big{)}$, vanishes thereby resulting in a
discontinuous drop of the bulk modulus. Conversely, in the compatible state
sufficient dilation will yield sudden contribution from perimeter tension. In
the following section, we apply a Landau-type energy analysis to understand
the trigger of shape zero-modes under dilation and compression.

##  IV Landau energy expansion

To understand how compression, or dilation, may trigger, or lift, shape
degeneracy we treat $\theta$ as an order parameter for the onset of shape
degeneracy such that $\theta\neq 0$ when cells adjust their shape to
accommodate imposed strains, and $\theta=0$ when cells remain rigid (no cell
response).

Before minimization, we expand the energy in power of $\theta$ to quartic
order,

$\displaystyle E(\epsilon,\theta;s_{0},r)\approx
E_{0}+\frac{\alpha}{2}(\epsilon,s_{0},r)\theta^{2}+\frac{\beta(\epsilon,s_{0},r)}{4}\theta^{4}+\mathcal{O}(\theta^{6})$
(9)

Where

$\displaystyle\alpha$
$\displaystyle=\frac{3^{\frac{5}{4}}}{2\sqrt{2}}r\ell(1+\epsilon)\bigg{(}s_{0}^{*}(6)\ell(1+\epsilon)-s_{0}\bigg{)},$
(10) $\displaystyle\beta$
$\displaystyle=\frac{3}{128}3^{\frac{1}{4}}r\ell(1+\epsilon)\bigg{(}13\sqrt{2}s_{0}-28\times
3^{\frac{1}{4}}(1+\epsilon)\ell\bigg{)}$ (11)

For the quartic approximation energy minimization yields
$\theta_{\text{min}}=0$ for $\alpha>0$, or
$\theta_{\text{min}}=\pm\sqrt{\frac{-\alpha}{\beta}}$ for $\alpha<0$. The
Landau expansion highlights the role of strain $\epsilon$ as tuning parameter
between the cell responding with $\theta_{\text{min}}=0$ or by spontaneously
tilting via a shear of $\theta_{\text{min}}=\pm\sqrt{\frac{-\alpha}{\beta}}$.
The form of $\alpha$ reflects an asymmetric response between compression
versus dilation. From $\alpha$ we can extract the critical strain,
$\epsilon_{*}$, which controls the onset/lifting of shape degeneracy.

$\displaystyle\epsilon_{*}=\frac{1}{\ell}\frac{s_{0}}{s_{0}^{*}(6)}-1$ (12)

The critical strain vanishes at the critical shape index and coincides with
the the failure of linear elasticity for any applied strain [3]. In the
compatible regime $\ell=1$ because target area is always achieved and
$\epsilon_{*}$ is independent of rigidity ratio $r$. Whereas in the
incompatible $\epsilon_{*}$ is dependent on $r$ through $\ell$.

We input the cell response via $\theta_{\text{min}}$ into the energy and
expand in powers of strain $\epsilon$.

$\displaystyle E(\epsilon;s_{0},r)=$
$\displaystyle\min_{\theta}E(\epsilon,\theta;s_{0},r)$ (13)
$\displaystyle\approx$
$\displaystyle\min_{\theta}E_{0}+\frac{\alpha}{2}(\epsilon,s_{0},r)\theta^{2}+\frac{\beta(\epsilon,s_{0},r)}{4}\theta^{4}+\mathcal{O}(\theta^{6})$
(14) $\displaystyle=$ $\displaystyle
E_{0}+\frac{1}{2}\left(\frac{\partial^{2}E(\theta_{\text{min}})}{\partial\epsilon^{2}}\right)\epsilon^{2}+\mathcal{O}(\epsilon^{3})$
(15)

In the final line the harmonic coefficient contains contributions from
$\theta$ which reduce the overall response of the tissue. If we did not
minimize over $\theta$ before expanding in $\epsilon$, the resultant deformed
energy would not reflect the self-shear response due to cell shape changes.

Figure 2: Plot A shows the mean field non-linear response bulk modulus versus
compression in the incompatible regime. At a critical strain $\epsilon_{*}$,
there is sudden discontinuous softening. The critical strain depends on
rigidity ratio, and we plot curves for $r=1,10,100$. In the compatible regime,
plot B shows the mean field bulk modulus hardening at a critical compression.
Note that $\epsilon_{*}$ in the incompatible regime is dependent on $r$ and
thereby sensitive to the balance between perimeter and areal elasticity.
Whereas in the compatible regime, the critical strain is only a function of
shape index. Plot C plot the difference in the linear response between
dilation and compression, and shows the continuous extension of the response
asymmetry extends away from the critical shape index. Plots A,B, and C
corresponds to hexagons. Plot D shows the effective critical shape index for a
random tiling of N cells on a sphere of radius $L=\sqrt{\frac{N}{4\pi}}$. The
mean field prediction is for pentagons , whereas the simulation is for a
disordered VM taken from [10].

A summary of the consequences of $\theta$ on the response are as follows: In
the incompatible regime, the $\theta_{\text{min}}=0$ solution corresponds to
the linear response in the solid state (see Fig. 3), whereas
$\theta_{\text{min}}=\pm\sqrt{\frac{-\alpha}{\beta}}$ corresponds to the
softer renormalized nonlinear response at critical compression strain. On the
other hand, in the compatible regime, the linear response is always given by
$\theta_{\text{min}}=\pm\sqrt{\frac{-\alpha}{\beta}}$ which allows perimeter
tension to vanish. The hardening under finite dilation occurs at a critical
dilative strain $\epsilon_{*}$ and corresponds to a switch from
$\theta_{\text{min}}=\pm\sqrt{\frac{-\alpha}{\beta}}$ to
$\theta_{\text{min}}=0$, resulting in a higher response. This hardening
phenomena is due to the cell’s inability access degenerate ground states to
accommodate large dilation.

### Calculating new critical shape index

So far our mean-field model has predicted how compression/dilation controls
the onset/lifting of shape degeneracy. If one defines the rigidity of a solid
by the availability of cell level zero modes, than the mean field treatment
suggests that dilation and compression shift the rigidity transition of the
VM.

The precise shift of the critical shape index can be found by solving
$\alpha=0$ for $s_{0}$ and setting the strain, $(1+\epsilon)$, as a material
parameter for local cell dilation/compression. Then solving for $s_{0}$
defines a ”new” critical shape index as set by dilation/compression. The
effective critical point $s_{0}^{\text{critical}}$ is

$\displaystyle\alpha$
$\displaystyle=\bigg{(}s_{0}^{*}(n)\ell(1+\epsilon)-s_{0}^{\text{critical}}\bigg{)}=0$
$\displaystyle\implies s_{0}^{\text{critical}}=s_{0}^{*}(n)\ell(1+\epsilon)$
(16)

Where we used the modified version of Eq. 10 for n-gons (see appendix C). Of
course this is not a complete shift of the rigidity transition since even if
perimeter tension vanishes there will always be areal tension. Hence the
system remains energetically frustrated, albeit less so, but soft.
Furthermore, $\epsilon$ needs to be determined by some cell growth/compression
process.

If cells are modeled on a 2D surface VM consisting of a 2D tiling embedded in
3D then cells may accommodate in-plane dilation/compression by buckling and
achieve true geometric compatibility in a curved geometry. For example,
geometrically frustrated (rigid) cells under compression may buckle to achieve
target area and perimeter by utilizing additional d.o.f. associated with
curvature, and thereby become compatible. By the same token, geometrically
compatible (floppy) cells experiencing dilation/growth will eventually
rigidity and become frustrated but may buckle to find a more
optimally/energetically lower frustrated state. Both case correspond to cells
utilizing curvature-related d.o.f to satisfy area and perimeter constraints
and thereby fine tuning their rigidity.

In support of this perceptive, simulations have reported that the rigidity
transition - as signalled by the critical target shape index - is indeed
sensitive to curvature. For example, a VM constrained on a sphere [10] has a
downward shifted $s_{0}^{*}$ dependent on curvature magnitude and relative
cell size.

In the rest of this article, we calculate the effective dilation/compression
set by curvature by calculating cell area on a curved surface relative to its
flat counterpart. Utilizing Eq.16, we predict the target shape index on curved
2D tissues.

##  V Rigidity transition in the presence of curvature.

We extend the mean field calculation by pertubatively calculating expressions
for cell area on surfaces of constant curvature in powers of
$K_{G}R_{\text{cell}}^{2}$, where $K_{G}$ is the Gaussian curvature, and
$R_{\text{cell}}$ the cell radius. Polygons and discs of fixed radius have
different total area depending on the surface on which they are inscribed. The
effective strain stems from the mismatch in area of curved cells from their
flat analog.

All geometric information of a surface $M$ is encoded in the metric tensor
$\mathbf{g}$ and we denote a surface generically as $(\mathbf{g},M)$. For a
general shape/cell on a surface, $D\subset M$, the area is defined by

$\displaystyle A$ $\displaystyle=\int_{D}\sqrt{\det g}d^{2}x$ (17)

Note that the determinant of the metric serves as a weight which accounts for
the local compression/dilation of the distances between points. Unlike the
planar setting of the mean field model where $g_{ij}=\delta_{ij}$, the metric
on curved surfaces - even uniformly curved - is not homogeneous but a function
of space. Nonetheless, the metric always admits a local pertubative expansion
in normal coordinates, $(x_{0},x^{i})$ which are defined by the condition that
geodesics can be locally parameterized as straight lines, i.e.
$\gamma(\lambda)=(x^{1}\lambda,x^{2}\lambda)$. In these coordinates, a
pertubative expansion of the metric in powers of curvature yields

$\displaystyle\det(g)=1-\frac{1}{3}\mathcal{R}|x|^{2}+\mathcal{O}(|x|^{3})$
(18)

Where $\mathcal{R}$ is the scalar curvature which is twice the Gaussian
curvature $K_{G}=2\mathcal{R}$. The expansion reflects how variations of the
metric are tied to curvature, and is locally approximated as flat with
corrections due to curvature. An outline of the derivation for Eq. 18 is given
in appendix A.

In our calculation we restrict attention to surfaces of uniform curvature by
fixing $\mathcal{R}$, and hence we only consider uniformly flat, sphereical,
and saddle-like surfaces. Of course, real curved biological tissues are not
not uniform and might have a boundary. Our approximation is controlled by the
dimensionless number set by the radius $R_{\text{cell}}$ of the cell over the
radius of curvature $L_{K}$

$\displaystyle C\equiv\frac{R_{\text{cell}}}{|L_{K}|}$ (19)

where $K_{G}\equiv\pm\frac{1}{L_{K}^{2}}$. Our mean field result will hold
best for tissues with moderate curvature or relativity small cells.

## Computing area

Using the metric expansion the area is expanded to quadratic order,

$\displaystyle A$ $\displaystyle=\int_{D}\sqrt{\det g}d^{2}x$
$\displaystyle=\int_{D}d^{2}x-\frac{\mathcal{R}}{3}\int_{D}|x|^{2}d^{2}x+\text{h.o.t.}$
(20)

The first term yields the flat case. We generalize our calculation for n-sided
polygons for easy comparison of various tilings.

To parameterize the polygonal n-sided cell $D$ we decompose it into 2n
triangles about the centroid as illustrated in figure 3. Details of the
calculation are given in appendix B. To quadratic order the area is

$\displaystyle
A=\bar{A}\left(1-\frac{\mathcal{R}}{6}R_{\text{cell}}^{2}f(n)+\mathcal{O}(R_{\text{cell}}^{4})\right)$
(21)

Where
$f(n)\equiv\cos^{2}\left(\frac{\pi}{n}\right)\left(\frac{2}{3}+\frac{1}{3}\sec^{2}\left(\frac{\pi}{n}\right)\right)$,
and
$\bar{A}=nR_{\text{cell}}^{2}\cos\left(\frac{\pi}{n}\right)^{2}\tan\left(\frac{\pi}{n}\right)$
is the flat area. In the limit of either very small cell size or very small
curvature, we get the typical planar area. Morever, in the limit
$n\rightarrow\infty$ the first correction yields
$\frac{\mathcal{R}}{6}R_{\text{cell}}^{2}\pi$ which reproduces the classical
result of Bertand-Diguet-Puiseux on the area comparison of a 2D geodesic ball
of radius $R_{\text{cell}}$ to its flat counterpart [12].

## Rigidity transition shift

From Eq.21 (and using Eq. 3) we write down the induced dilation/compression
strain set by curvature

$\displaystyle\epsilon_{\mathcal{R}}(n)$
$\displaystyle=-1+\sqrt{1-\frac{\mathcal{R}}{6}R_{\text{cell}}^{2}f(n)}$
$\displaystyle=-1+\sqrt{1-\cos^{2}\left(\frac{\pi}{n}\right)\left(\frac{2}{3}+\frac{1}{3}\sec^{2}\left(\frac{\pi}{n}\right)\right)\frac{\mathcal{R}}{12}R^{2}_{\text{cell}}}$
(22)

The effective strain depends on the number of edges due to the discrete
rotational symmetry of polygons: the furthest edges are weighted differently
than the midpoint of edges. This vanishes in the $n=\infty$ limit as the cell
is approximates a disc and regains full rotational symmetry.

The predicted shift in the critical shape index is

$\displaystyle s_{0}^{\text{critical}}(n,\mathcal{R})$
$\displaystyle=s_{0}^{*}(n)\ell(n)(1+\epsilon_{\mathcal{R}}(n))$ (23)
$\displaystyle\approx s_{0}^{*}(n)(1+\epsilon_{\mathcal{R}}(n))$

Where we’ve set $\ell\approx 1$, which restricts our prediction near the
planar critical point or large rigidity ratio $r>>1$.

## Comparison with simulation

In the work by Sussman [10] a disordered vertex model of $N$ cells on a
uniform sphere of radius $L=\sqrt{\frac{N}{4\pi}}$ was simulated for various
$N$. At the onset of rigidity, the observable shape index per cell
$\frac{P}{\sqrt{A}}$ was extracted and used to compute a numerical probability
distribution. The peak of this distribution was used to define the critical
shape index on the sphere. For large systems, the critical point was found to
be sharply peaked at around $s_{0}^{*}\sim 3.812$, which corresponds to the
the critical shape index for pentagons[13].

To compare with [10] we re-cast $\mathcal{R}R_{\text{cell}}^{2}$ in terms of
$N$. The radius of curvature for a sphere gives the Gaussian curvature
$K_{G}\equiv\frac{1}{L^{2}}$, which is 1/2 the scalar curvature.

$\displaystyle s_{0}^{\text{critical}}(n=5,K_{G})\approx
3.812\left(1-\frac{(7+\sqrt{5})}{72}\frac{2R_{\text{cell}}^{2}}{L^{2}}\right)^{1/2}$
(24)

Sussman works on a sphere of raduis $L=\sqrt{\frac{N}{4\pi}}$ and average cell
area is set to unity, i.e. $A_{\text{cell}}\equiv\frac{N}{4\pi L^{2}}=1$. We
neglect cell packing since our mean field calculation is for a single cell,
and so we take $R_{\text{cell}}$ as given for a pentagon. Thus
$5R_{cell}^{2}\cos\left(\frac{\pi}{5}\right)^{2}\tan\left(\frac{\pi}{5}\right)=1\implies$
$R_{\text{cell}}^{2}\approx\frac{1}{2.377}$.

Therefore the relative ratio of cell size to radius of curvature goes as
$\frac{R_{\text{cell}}^{2}}{L^{2}}\approx\frac{4\pi}{2.377}\frac{1}{N}$,
yielding

$\displaystyle s_{0}^{\text{critical}}(n=5,K_{G})\approx
3.812\left(1-\frac{(7+\sqrt{5})}{36}\frac{4\pi}{2.377}\frac{1}{N}\right)^{1/2}$
(25)

Comparison to simulation data is shown in Fig. 2 Plot D.

Besides expanding to higher order the calculation can be improved by computing
the ground state characteristic cell size $\ell_{0}$ for curved vertex models,
but this calculation is beyond our mean field approach. Additionally, [10]
reports that the shape index distribution broadens for larger
$R^{2}_{\text{cell}}K_{G}$, reflecting a greater diversity of polygons at the
rigidity transition than the flat counterpart. Taking into account this
greater diversity could help refine the curvature correction in Eq. 25. In
particular, for relative large curvature other polygonal shapes besides the
pentagon could more relevant for disordered systems.

##  VI Discussion

Utilizing a mean field model we show the asymmetry of the linear response in
the vertex model under dilation and compression extends away from the critical
shape index for finite strains. The asymmetry reflects how an initially rigid
tissue may be sufficiently compressed to induce shape degeneracy and thereby
relax perimeter tension, yielding a softer bulk modulus. Furthermore,
sufficient dilation applied to a compatible (floppy) cell lifts shape
degeneracy, yielding an uptick in the bulk modulus. Thus applied dilation and
compression shift the rigidity of the VM in 2D.

Using this insight, we extend our mean field theory to calculate the effective
dilation/compression engendered by intrinsic curvature and predict the precise
shift in the rigidity transition by calculating the effective critical shape
index. We compare our result to simulation done by [10] and find good
agreement.

Our mean field prediction provides a metric which can be applied to studying
the rigidity transition in curved biological tissues in a manner analogous as
the planar shape index.

##  VII Appendices

###  VII.1 Details about metric expansion

The purpose of this appendix is to give a brief explanation of the series
expansion of the metric in terms of curvature. A more complete and rigorous
treatment may be found in many textbooks on Riemannain geometry such as [14,
15], and in particular we follow the classic treatise [16].

In general the metric is a second order tensor whose components are spatially
dependent function of the surface. It governs all geometric data in that the
distance between any two points is given by the line element

$\displaystyle ds^{2}=g_{ij}(x)dx^{i}dx^{j}.$ (26)

About a given point $x_{0}\in M$ in some neighborhood, the components of the
metric are approximated as constants to 1st order. One may diagonalize this
linear approximation such that the metric at $x_{0}$ is given by
$\delta_{ij}$. However to 2nd order the metric’s components need not be also
constants. In fact, if there exist coordinates such that the metric’s Taylor
series is constant up to 2nd order then the metric is totally flat in the
neighborhood.

Normal coordinates about a point $x_{0}$ are defined locally where geodesics
may be parametrized in local coordinates $x^{i}$ by $\lambda$ such that
$\gamma^{i}(\lambda)=x^{i}\lambda$, where $\gamma(0)\equiv x_{0}$. The
Christoffel symbols are extracted from the geodesic equation

$\displaystyle
0=\frac{d^{2}\gamma^{i}}{d\lambda^{2}}+\Gamma_{k\ell}^{i}\frac{d\gamma^{k}}{d\lambda}\frac{d\gamma^{\ell}}{d\lambda}$
(27)

Utilizing normal coordinates, Eq. 29 immediately implies
$\Gamma_{k\ell}^{i}(x_{0})=0$. Differentiation and index manipulation also
yields the differential constraint equation.

$\displaystyle\partial_{j}\Gamma_{k\ell}^{i}(x_{0})+\partial_{k}\Gamma_{\ell
j}^{i}(x_{0})+\partial_{\ell}\Gamma_{jk}^{i}(x_{0})=0$ (28)

The Riemann curvature tensor is defined as

$\displaystyle R^{i}_{jkl}$
$\displaystyle=\partial_{k}\Gamma^{i}_{jl}-\partial_{l}\Gamma^{i}_{jk}+\Gamma^{i}_{pk}\Gamma^{p}_{kl}+\Gamma^{i}_{pl}\Gamma^{p}_{kj}$
(29)

From the differential constraint and the definition of $R^{i}_{jkl}$, one can
show

$\displaystyle\partial_{l}\Gamma_{ij}^{k}=-\frac{1}{3}\left(R_{ijl}^{k}+R_{jil}^{k}\right)$
(30)

Symmetry of the metric implies the covariant derivative of the metric
vanishes,i.e. $\nabla\mathbf{g}=0$
$\implies\partial_{k}g_{ij}-\Gamma^{p}_{jk}g_{ip}-\Gamma^{p}_{ik}g_{jp}=0$.
The second derivative of the metric in thes coordinates is

$\displaystyle\partial_{kl}^{2}g_{ij}=-\frac{1}{3}\left(R_{klij}+R_{jlik}\right),$
(31)

The Taylor expansion of the metric in normal coordinates yields

$\displaystyle g_{ij}$
$\displaystyle=\delta_{ij}-\frac{1}{3}R_{ijk\ell}x^{k}x^{\ell}+\mathcal{O}(|x|^{2})$
(32)

Higher order terms can be generated iteratively by calculating higher order
differential constraint equations from Eqs.28 and $\nabla\mathbf{g}=0$.

For 2D surfaces the Riemann curvature tensor only has a single d.o.f. and
admits the representation [15]

$\displaystyle R_{iklj}=\mathcal{R}(g_{ik}g_{lj}-g_{ij}g_{kl})$ (33)

Where $\mathcal{R}=2K_{G}$ is the Ricci scalar curvature. From this the Ricci
tensor follows $R_{ij}\equiv g^{kl}R_{iklj}=\mathcal{R}g_{ij}$. Using the
expansion of the metric, we have to lowest order

$\displaystyle R_{iklj}$
$\displaystyle=\mathcal{R}(\delta_{ik}\delta_{lj}-\delta_{ij}\delta_{kl})+\mathcal{O}(|x|^{2})$
(34) $\displaystyle R_{ij}$
$\displaystyle=\mathcal{R}\delta_{ij}+\mathcal{O}(|x|^{2})$ (35)

These expressions reflect that locally any surface looks either flat
$(\mathcal{R}=0)$, spherical $(\mathcal{R}>0)$, or saddle-like
$(\mathcal{R}<0)$.

To lowest order the metric expansion about $p$ becomes

$\displaystyle g_{ij}=$
$\displaystyle\delta_{ij}-\frac{1}{3}\mathcal{R}(\delta_{ik}\delta_{lj}-\delta_{ij}\delta_{kl})x^{k}x^{l}+\mathcal{O}(|x|^{3})$
(36)

and determinant yields

$\displaystyle\det(g)=1-\frac{1}{3}\mathcal{R}|x|^{2}+\mathcal{O}(|x|^{3})$
(38)

which shows explicitly how curvature induces local compression or dilation.
Higher order terms contain gradients and higher order invariants of
$R_{ijkl}$, and therefore are completely determined by $K$. Thus it follows
that if the quadratic contribution vanishes, then the metric is totally flat
in the neighborhood.

###  VII.2 Pertubative polygon area expansion

To explicitly parameterize the polygonal cell $D$, we will consider a regular
n-gon and decompose it into 2n-triangles about its centroid as pictured in
figure 3.

Working in terms of polar coordinates, this yields for the first term

$\displaystyle\int_{D}d^{2}x$
$\displaystyle=2n\int_{0}^{\frac{\pi}{n}}\int_{0}^{R_{\text{cell}}\cos(\frac{\pi}{n})\sec\theta}d\theta
rdr$
$\displaystyle=nR_{\text{cell}}^{2}\cos\left(\frac{\pi}{n}\right)^{2}\tan\left(\frac{\pi}{n}\right)$
(39)

In the limit of $n\rightarrow\infty$ we get $\pi R_{\text{cell}}^{2}$, as
expected for circles. Using the same coordinate system, we compute the first
correction due to curvature

$\displaystyle\frac{\mathcal{R}}{3}\int_{D}|x|^{2}d^{2}x$
$\displaystyle=\frac{\mathcal{R}}{3}2n\int_{0}^{\frac{\pi}{n}}\int_{0}^{R_{\text{cell}}\cos(\frac{\pi}{n})\sec\theta}d\theta
r^{3}dr$
$\displaystyle=\frac{\mathcal{R}}{6}R_{\text{cell}}^{4}\cos^{4}\left(\frac{\pi}{n}\right)n\left(\frac{2}{3}+\frac{1}{3}\sec^{2}\left(\frac{\pi}{n}\right)\right)\tan\left(\frac{\pi}{n}\right)$
(40)

###  VII.3 Mean field vertex model

The mean field model is defined by the area and perimeter of a single cell,
which is parameterized by n-edges $\nu_{\alpha}$ given by

$\displaystyle\vec{\nu}_{\alpha}\equiv\ell_{0}\left(\cos\left(\frac{2\pi\alpha}{n}\right),\sin\left(\frac{2\pi\alpha}{n}\right)\right)$
(41)

Where $\ell_{0}$ the characteristic cell edge length. The perimeter is the sum
of each edge length

$\displaystyle P=$
$\displaystyle\sum_{\alpha}^{n}\sqrt{\vec{\nu}_{\alpha}\cdot\vec{\nu}_{\alpha}}$
(42)

Under an affine transformation, denoted as the matrix $\mathbf{F}$, the
perimeter transformations as

$\displaystyle P=$
$\displaystyle\sum_{\alpha}^{n}\sqrt{(\mathbf{F}\vec{\nu}_{\alpha})\cdot(\mathbf{F}\vec{\nu}_{\alpha}})$
(43)

The area can be calculated by the cross product

$\displaystyle A=\int_{D}dx^{2}=n|\vec{a}\times\vec{b}|$ (44)

where $\vec{a}$ and $\vec{b}$ are defined in Fig.3.

The polygon will be changed under some affine transformation $\mathbf{F}$. The
area term is striaigtforeward since cross product transforms as
$|(\mathbf{F}\vec{a})\times(\mathbf{F}\vec{b})|=\det(\mathbf{F})|\vec{a}\times\vec{b}|$.
Thus under any affine transformation, the area term takes the simple form

$\displaystyle A=$ $\displaystyle\det(\mathbf{F})n|\vec{a}\times\vec{b}|$ (45)
$\displaystyle=$
$\displaystyle\det(\mathbf{F})\frac{n}{4}\ell_{0}^{2}\cot\left(\frac{\pi}{n}\right)$
(46)

Figure 3: To make our calculation easily applicable to squares, pentagons, and
hexagons we work with this figure to switch between the different polygons and
different choices of encoding cell shape either via edge length $\ell_{0}$ or
cell radius $R_{\text{cell}}$.

The energy per cell is cast as

$\displaystyle
E=\frac{\kappa_{A}}{2}\left(\frac{n}{4}\ell_{0}^{2}\cot\left(\frac{\pi}{n}\right)\det(\mathbf{F})-A_{0}\right)^{2}$
$\displaystyle+\frac{\kappa_{P}}{2}\left(\sum_{\alpha}^{n}\sqrt{(\mathbf{F}\vec{\nu}_{\alpha})\cdot(\mathbf{F}\vec{\nu}_{\alpha}})-P_{0}\right)^{2}$
(47)

To non-dimensionalize we define reference lengths $\ell_{A}$ and $\ell_{P}$
such that

$\displaystyle A_{0}$
$\displaystyle=\frac{n}{4}\ell_{A}^{2}\cot\left(\frac{\pi}{n}\right)$ (48)
$\displaystyle P_{0}$ $\displaystyle=n\ell_{P}$ (49)

And rescale energy by $\kappa_{A}A_{0}^{2}$, yielding

$\displaystyle
E=\frac{1}{2}\left(\ell^{2}\det(\mathbf{F})-1\right)^{2}+\frac{r}{2}\left(\ell\sum_{\alpha}^{n}\sqrt{(\mathbf{F}\vec{\nu}_{\alpha})\cdot(\mathbf{F}\vec{\nu}_{\alpha}})-s_{0}\right)^{2}$
(50)

Where $r\equiv\frac{\kappa_{P}}{\kappa_{A}A_{0}}$,
$s_{0}\equiv\frac{P_{0}}{\sqrt{A_{0}}}$ is the target shape index, and
$\ell\equiv\frac{\ell_{0}}{\sqrt{A_{0}}}$ is the re-scaled characteristic cell
edge length.

In the incompatible state, the ground state corresponds to a regular polygon
with $\ell_{0}$ defined to minimize the energy. This involves solving the
following cubic equation.

$\displaystyle\frac{\partial E}{\partial\ell}=0$ (51)

The relevant solution obeys the inequality $\ell\leq 1$ for all $s_{0}\leq
s_{0}^{*}$. In the compatible state energy minimization yields that
$\ell\equiv 1$ for choices of $r$ and $s_{0}\geq s_{0}^{*}$.

#### Inputting deformations

The linear transformation $\mathbf{F}$ encodes both applied deformations and
cell response. We assume that all shape distortions of the cell can be
captured in by linear affine transformation. Although the non-affine
contribution is important in calculating exactly even the linear response of
the vertex model as shown by the authors [11].

The applied deformation is set by a transformation rule. In this article we
focus on compression and dilation, which yields the matrix

$\displaystyle\mathbf{F}_{\epsilon}=\begin{pmatrix}1+\epsilon&0\\\
0&1+\epsilon\end{pmatrix}$ (52)

On top of this, we also will allow the cell to adjust its perimeter without
changing the imposed re-scaled area. This imposes the constraint

$\displaystyle\det(\mathbf{F}^{\text{cell}})=1$ (53)

Which only fixes a single degree of freedom, leaving in principle three
components of $\mathbf{F}_{\text{cell}}$ free. For simplicity, we only
consider the the cell’s response by tilting through a simple shear
transformation

$\displaystyle\mathbf{F}^{\text{cell}}_{\theta}=\begin{pmatrix}1&\tan(\theta)\\\
0&1\end{pmatrix}$ (54)

Note that our mean field model is invariant under rotations since area and
perimeter are rotationally invariant objects. Additionally, if we include
other affine transformation, all their respectively matrices $\mathbf{F}$ will
commute with one another.

We set the overall net deformation gradient in the mean field model as

$\displaystyle\mathbf{F}=\mathbf{F}^{\text{cell}}_{\theta}\cdot\mathbf{F}_{\epsilon}=\begin{pmatrix}1&\tan(\theta)\\\
0&1\end{pmatrix}\begin{pmatrix}1+\epsilon&0\\\ 0&1+\epsilon\end{pmatrix}$ (55)

With these two shape degrees of freedom, $\epsilon$ and $\theta$, we study the
non-linear response to finite compression/dilation and extract our rigidity
transition shift due to curvature.

## References

  * Park  _et al._ [2015] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T. Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H. Kim, B. Gweon, _et al._ , Unjamming and cell shape in the asthmatic airway epithelium, Nature materials 14, 1040 (2015).
  * Farhadifar [2009] R. Farhadifar, Dynamics of cell packing and polar order in developing epithelia, PhD thesis (2009).
  * Hernandez  _et al._ [2021] A. Hernandez, M. F. Staddon, M. J. Bowick, M. C. Marchetti, and M. Moshe, Geometric rigidity and anomalous elasticity of cellular tissue vertex model, arXiv preprint arXiv:2109.10407 (2021).
  * Huang  _et al._ [2021] J. Huang, J. O. Cochran, S. M. Fielding, M. C. Marchetti, and D. Bi, Shear-driven solidification and nonlinear elasticity in epithelial tissues, arXiv preprint arXiv:2109.10374 (2021).
  * Bi  _et al._ [2015] D. Bi, J. Lopez, J. M. Schwarz, and M. L. Manning, A density-independent rigidity transition in biological tissues, Nature Physics 11, 1074 (2015).
  * Osserman  _et al._ [1978] R. Osserman _et al._ , The isoperimetric inequality, Bulletin of the American Mathematical Society 84, 1182 (1978).
  * Farhadifar  _et al._ [2007] R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing, Current Biology 17, 2095 (2007).
  * Staple  _et al._ [2010] D. B. Staple, R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher, Mechanics and remodelling of cell packings in epithelia, The European Physical Journal E 33, 117 (2010).
  * Murisic  _et al._ [2015] N. Murisic, V. Hakim, I. G. Kevrekidis, S. Y. Shvartsman, and B. Audoly, From discrete to continuum models of three-dimensional deformations in epithelial sheets, Biophysical journal 109, 154 (2015).
  * Sussman [2020] D. M. Sussman, Interplay of curvature and rigidity in shape-based models of confluent tissue, Physical Review Research 2, 023417 (2020).
  * Staddon  _et al._ [2022] Staddon, Hernandez, Bowick, Moshe, and Marchetti, Elastic moduli of the vertex model, manuscript in preparation (2022).
  * Gray and Vanhecke [1979] A. Gray and L. Vanhecke, Riemannian geometry as determined by the volumes of small geodesic balls, Acta Mathematica 142, 157 (1979).
  * Bi  _et al._ [2016] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, Motility-driven glass and jamming transitions in biological tissues, Physical Review X 6, 021011 (2016).
  * Gray [2003] A. Gray,  _Tubes_ , Vol. 221 (Springer Science & Business Media, 2003).
  * Willmore [1993] T. Willmore,  _Riemannian Geometry_ (Oxford Science Publications, 1993).
  * Veblen and Whitehead [1932] O. Veblen and J. Whitehead,  _Foundations of Differential Geometry_ (Cambridge Press, 1932).


# Carrier Frequency Offset Estimation for OCDM with Null Subchirps

Sidong Guo, , Yiyin Wang, and Xiaoli Ma

###### Abstract

In this paper, we investigate the carrier frequency offset (CFO)
identifiability problem in orthogonal chirp division multiplexing (OCDM)
systems. We propose a transmission scheme by inserting consecutive null
subchirps. A CFO estimator is accordingly developed to achieve a full
acquisition range. We further demonstrate that the proposed transmission
scheme not only help to resolve CFO identifiability issues but also enable
multipath diversity for OCDM systems. Simulation results corroborate our
theoretical findings.

###### Index Terms:

CFO, OCDM, multipath diversity, identifiability

##  I Introduction

Orthogonal chirp division multiplexing (OCDM) is a multi-carrier modulation
scheme that uses orthogonal chirp waveforms to carry information symbols. OCDM
originates from paraxial optics and is first applied to wireless
communications [1] [2]. Compared to the currently dominant modulation scheme,
i.e., orthogonal frequency division multiplexing (OFDM), OCDM is more robust
against frequency-selective channels, enjoys lower bit error rate (BER) under
various equalization methods, and performs better against burst interferences.
Gaining more attractions recently, OCDM also sees increased interests in the
research of multi-user OCDM, channel and CFO estimation [3, 4, 5, 6, 7].

In an OCDM system, the inverse discrete Fresnel transform (IDFnT) is used as
the transform kernel to spread each information symbol over the entire
bandwidth and the whole symbol block. Thus, OCDM has the double-spreading
feature. A transform basis being non-sinusoidal may enjoy certain advantages
against CFO. For instance, simulation results in [8] indicate that OCDM with a
minimum mean square error (MMSE) equalizer or decision feedback equalizer
(DFE) has a much better performance in terms of BER than OFDM over frequency-
selective channels with CFO at medium signal-to-noise ratios (SNRs). It is
also observed in [9] that chirp-based OFDM systems outperforms traditional
OFDM with the presence of CFO.

The double-spreading property of chirp basis suggests that frequency-selective
channels may have different impacts on receiver’s ability to correctly
identify CFO in OCDM compared to OFDM. The identifiability issue is well
documented in the case of OFDM in [10]. Owing to the sinusoidal kernel, OFDM
converts a frequency-selective channel into a set of frequency-flat channels.
However, channel nulls in the frequency domain can create multiple
indistinguishable candidate CFOs in OFDM systems. Thus, it was proposed that
judicially placed null subcarriers can restore CFO identifiability in OFDM,
where channel nulls are sufficiently covered by null subcarriers [10, 11]. The
optimality and identifiability of data-aided CFO estimation approaches for
OFDM systems with training sequences are explored in [12, 13].

On the other hand, some preliminary efforts are directed towards handling CFO
in OCDM. In [5], an excessive cyclic prefix (CP) is used to cross correlate
with its counterpart in an OCDM block to estimate the CFO. In [6], a preamble-
type synchronization and channel estimation method is applied to an OCDM
system. Although both of them [5, 6] show decent accuracy, they require
redundant CPs or symbols. They also face identifiability issues at a high
normalized CFO. For instance, the CP-based approach [5] has ambiguity in the
CFO estimation, when the CFO exceeds half of the subchirp frequency spacing.

This paper resolves the CFO identifiability issue in OCDM with the help of the
proposed transmission scheme, which inserts consecutive null subchirps. A CFO
estimator is designed accordingly to achieve a full acquisition range of CFO.
We further prove that the proposed transmission scheme not only help to
address CFO identifiability issue, but also enable multipath diversity for
OCDM systems. Simulation results verify the uniqueness of the proposed CFO
estimator, and indicate that the BER performances of both the linear
equalizers (LEs) and maximum likelihood estimator (MLE) benefit from the
enabled multipath diversity.

Notations: Upper (lower) case bold letters are used to represent matrices
(vectors). The notation $[\mathbf{A}]_{m,n}$ indicates the entry at the $m$th
row and $n$th column of the matrix $\mathbf{A}$, whereas $[\mathbf{A}]_{m,:}$
and $[\mathbf{A}]_{:,n}$ denotes the $m$th row and $n$th column of the matrix
$\mathbf{A}$, respectively. The identity matrix of size $P\times P$ is denoted
by $\mathbf{I}_{P}$. The zero matrix of size $M\times N$ is denoted by
$\mathbf{0}_{M\times N}$. The notations $(\cdot)^{T}$,
$(\cdot)^{\mathcal{H}}$, $(\cdot)^{-1}$, $(\cdot)^{\dagger}$, and
$(\cdot)^{\ast}$ represent transpose, conjugate transpose, inverse, Moore-
Penrose inverse, and element-wise conjugate, respectively. We use $E[\cdot]$
to denote the expectation with respect to (w.r.t.) all random variables within
the brackets.

##  II System Model

In this section, we derive a baseband equivalent OCDM system model under
frequency-selective channels with CFO, where the time-invariant frequency-
selective channels are assumed with channel impulse responses (CIRs) $h(l)$ of
order $L$, where $l=0,1,\dots,L$.

Consider the $i$th symbol block is composed of the $i$th data symbol block
$\mathbf{s}(i)$ and null symbols. The null symbols in the Fresnel domain
correspond to null subchirps in the time domain. The null symbols are added in
a consecutive way. Thus, the $i$th symbol block of length $N$ is assembled as
$\mathbf{T}_{\rm zp}\mathbf{s}(i)$, where $\mathbf{s}(i)$ is of length $K$ and
$\mathbf{T}_{\rm zp}=[\mathbf{I}_{K}\quad\mathbf{0}_{K\times(N-K)}]^{T}$. We
further assume that the covariance matrix of $\mathbf{s}(i)$ is given by
$\mathbf{R}_{ss}=E[\mathbf{s}(i)\mathbf{s}^{\mathcal{H}}(i)]=E_{s}\mathbf{I}_{K}$,
where $E_{s}$ is the symbol energy. The $i$th symbol block is modulated with
the IDFnT. The resulting time-domain symbol block can be written as
$\mathbf{x}(i)=\mathbf{\Phi}^{\mathcal{H}}\mathbf{T}_{\rm zp}\mathbf{s}(i)$,
where the DFnT matrix $\mathbf{\Phi}$ of size $N\times N$ is defined as:

$\displaystyle[\mathbf{\Phi}]_{m,n}=\frac{1}{\sqrt{N}}e^{-j\frac{\pi}{4}}\times$
$\displaystyle\begin{cases}e^{j\frac{\pi}{N}(m-n)^{2}},&N=0\,\,{\rm
mod}\,\,2,\\\ e^{j\frac{\pi}{N}(m+\frac{1}{2}-n)^{2}},&N=1\,\,{\rm
mod}\,\,2,\end{cases}$ $\displaystyle m,n=1,\dots,N.$ (1)

The matrix $\mathbf{\Phi}$ is circulant and unitary, i.e.,
$\mathbf{\Phi}^{\mathcal{H}}\mathbf{\Phi}=\mathbf{\Phi}\mathbf{\Phi}^{\mathcal{H}}=\mathbf{I}_{N}$.
Moreover, a cyclic prefix (CP) of length $L$ is inserted in the time domain
before $\mathbf{x}(i)$ to combat delay effects.

The transmitted OCDM signal goes through a time-invariant frequency-selective
channel. At the receiver side, we consider a normalized CFO $w_{o}=2\pi
f_{o}T_{s}$ in the full range $[-\pi,\pi)$, where $f_{o}$ is the CFO in Hz,
and $T_{s}$ is the sampling period. After synchronization and removal of the
CP, the received signal block $\mathbf{y}(i)$ of length $N$ is given by

$\mathbf{y}(i)=e^{jw_{o}(i(N+L)-N)}\mathbf{D}_{N}(w_{o})\mathbf{H}\mathbf{\Phi}^{\mathcal{H}}\mathbf{T}_{\rm
zp}\mathbf{s}(i)+\mathbf{n}(i),$ (2)

where $\mathbf{D}_{N}(w_{o})$ of size $N\times N$ is the diagonal CFO matrix
with the $n$th element in the diagonal given by $e^{jw_{o}(n-1)}$, the first
column of the circulant channel matrix $\mathbf{H}$ is given by
$\mathbf{h}=[h(0),\dots,h(L),\mathbf{0}_{N-L-1}^{T}]^{T}$, and $\mathbf{n}(i)$
is the independent and identically distributed (i.i.d.) zero-mean additive
white Gaussian noise (AWGN) vector with the covariance matrix
$\sigma^{2}\mathbf{I}_{N}$.

Once the CFO estimate $\hat{w}_{o}$ is obtained, we compensate the CFO effect
on the received signal block and achieve
$\mathbf{r}(i)=e^{-j\hat{w}_{o}(i(N+L)-N)}\mathbf{D}_{N}^{\mathcal{H}}(\hat{w}_{o})\mathbf{y}(i)$.
Consequently, LEs (denoted by $\bf G_{(\cdot)}$) or non-linear equalizers can
be applied to $\mathbf{r}(i)$. For example, the data symbols can be recovered
using LEs as $\hat{\mathbf{s}}(i)=\mathbf{G}_{(\cdot)}\mathbf{r}(i)$, where
the zero-forcing (ZF) or MMSE equalizers are given by

$\displaystyle\mathbf{G}_{\rm ZF}$ $\displaystyle=\mathbf{B}^{\dagger},$ (3)
$\displaystyle\mathbf{G}_{\rm MMSE}$
$\displaystyle=\mathbf{B}^{\mathcal{H}}\left(\frac{\sigma^{2}}{E_{s}}\mathbf{I}_{N}+\mathbf{B}\mathbf{B}^{\mathcal{H}}\right)^{-1},$
(4)

with $\mathbf{B}={\mathbf{H}}\mathbf{\Phi}^{\mathcal{H}}\mathbf{T}_{\rm zp}$.

##  III CFO Estimator and its Identifiability

In this section, we develop an estimator for the true CFO $w_{0}$ using the
left null space (LNS) of the covariance matrix
$\mathbf{R}_{yy}=E[\mathbf{y}(i)\mathbf{y}^{\mathcal{H}}(i)]$. We first
analyze the LNS of ${\mathbf{R}}_{yy}$. Moreover, a CFO estimator based on the
LNS is proposed, and its uniqueness is proved.

###  III-A The Proposed CFO Estimator

Based on 2, we derive the covariance matrix $\mathbf{R}_{yy}$ as follows,

$\mathbf{R}_{yy}=\mathbf{D}_{N}(w_{0})\mathbf{\Phi}^{\mathcal{H}}\mathbf{H}\mathbf{T}_{\rm
zp}\mathbf{R}_{ss}\mathbf{T}_{\rm
zp}^{\mathcal{H}}\mathbf{H}^{\mathcal{H}}\mathbf{\Phi}\mathbf{D}_{N}^{\mathcal{H}}(w_{0})\\!+\\!\sigma^{2}\mathbf{I}_{N},$
(5)

where (5) is achieved using the property of circulant matrices
$\mathbf{\mathbf{H}\Phi}^{\mathcal{H}}=\mathbf{\Phi}^{\mathcal{H}}\mathbf{H}$.
Furthermore, the noiseless part of the covariance matrix is defined as

$\bar{\mathbf{R}}_{yy}=\mathbf{\Phi}^{\mathcal{H}}\mathbf{H}\mathbf{T}_{\rm
zp}\mathbf{R}_{ss}\mathbf{T}_{\rm
zp}^{\mathcal{H}}\mathbf{H}^{\mathcal{H}}\mathbf{\Phi},$ (6)

where $\bar{\mathbf{R}}_{yy}$ is of size $N\times N$ and ${\rm
rank}(\mathbf{T}_{\rm zp})=K$. Thus, the rank of $\bar{\mathbf{R}}_{yy}$ is
upper-bounded by ${\rm rank}(\bar{\mathbf{R}}_{yy})\leq{\rm min}\left({\rm
rank}({\mathbf{H}}),{\rm rank}(\mathbf{T}_{\rm zp})\right)<N$, according to
Sylvester’s inequality. Hence, there exists the LNS of
$\bar{\mathbf{R}}_{yy}$.

The zero-padding (ZP) matrix $\mathbf{T}_{\rm zp}$ makes
$\mathbf{H}\mathbf{T}_{\rm zp}$ a Toeplitz matrix of size $N\times K$. Recall
that the multipath channel order is $L$. When $N-K>L$, the number of inserted
null subchirps is greater than $L$, and the last $N-K-L$ rows of
$\mathbf{H}\mathbf{T}_{\rm zp}$ only contain zeros. As a result, we arrive at

$\mathbf{\Phi}^{\mathcal{H}}\mathbf{H}\mathbf{T}_{\rm
zp}=\left[\mathbf{\Phi}^{\ast}_{\mathcal{S}}\,\,\mathbf{\Phi}^{\ast}_{\mathcal{N}}\right]\left[\begin{pmatrix}\bar{\mathbf{H}}\\\
\mathbf{0}_{(N-K-L)\times K}\end{pmatrix}\right],$ (7)

where
$\mathbf{\Phi}^{\ast}_{\mathcal{S}}=\left[\mathbf{\Phi}^{\mathcal{H}}\right]_{:,1:K+L}$
of size $N\times(K+L)$,
$\mathbf{\Phi}^{\ast}_{\mathcal{N}}=\left[\mathbf{\Phi}^{\mathcal{H}}\right]_{:,K+L+1:N}$
of size $N\times(N-K-L)$, and
$\bar{\mathbf{H}}=\left[\mathbf{H}\right]_{1:K+L,1:K}$ of size $(K+L)\times
K$. According to 7, we further prove that the null subchirps
$\mathbf{\Phi}_{\mathcal{N}}=\left[\mathbf{\Phi}^{T}\right]_{:,K+L+1:N}$ is in
the LNS of $\bar{\mathbf{R}}_{yy}$ as

$\displaystyle\mathbf{\Phi}^{T}_{\mathcal{N}}\left[\mathbf{\Phi}^{\ast}_{\mathcal{S}}\,\,\mathbf{\Phi}^{\ast}_{\mathcal{N}}\right]\left[\begin{pmatrix}\bar{\mathbf{H}}\\\
\mathbf{0}_{(N-K-L)\times K}\end{pmatrix}\right]$ $\displaystyle=$
$\displaystyle\left[\mathbf{0}_{(N-K-L)\times(K+L)}\,\,\mathbf{I}_{N-K-L}\right]\left[\begin{pmatrix}\bar{\mathbf{H}}\\\
\mathbf{0}_{(N-K-L)\times K}\end{pmatrix}\right]$ $\displaystyle=$
$\displaystyle\mathbf{0}_{(N-K-L)\times K}.$ (8)

Define a set of indices $\mathcal{N}_{\phi}=\\{K+L+1,\dots,N\\}$, with this
notation, the CFO $w_{0}$ can be estimated by using the LNS of
$\bar{\mathbf{R}}_{yy}$ as

$\hat{w}_{0}=\mathop{\arg\min}\limits_{w}J(w),$ (9)

where the cost function $J(w)$ is defined as

$J(w)=\sum_{k\in\mathcal{N}_{\phi}}\boldsymbol{\phi}_{k}^{T}\mathbf{D}_{N}^{-1}(w){\mathbf{R}}_{yy}\mathbf{D}_{N}(w)\boldsymbol{\phi}_{k}^{\ast},$
(10)

with $\boldsymbol{\phi}_{k}=\left[\mathbf{\Phi}^{T}\right]_{:,k}$ and
$\boldsymbol{\phi}^{\ast}_{k}=\left[\mathbf{\Phi}^{\mathcal{H}}\right]_{:,k}$.
Leveraging $\bar{\mathbf{R}}_{yy}$, the cost function 10 can be rewritten as

$\displaystyle J(w)=$
$\displaystyle\sum_{k\in\mathcal{N}_{\phi}}\boldsymbol{\phi}_{k}^{T}\mathbf{D}_{N}(w_{0}-w)\bar{\mathbf{R}}_{yy}\mathbf{D}_{N}(w-w_{0})\boldsymbol{\phi}_{k}^{\ast}$
$\displaystyle+\sigma^{2}(N-K-L),$ (11)

where the equivalent noise term $\sigma^{2}(N-K-L)$ is the result of
$\sigma^{2}\sum_{k\in\mathcal{N}_{\phi}}\boldsymbol{\phi}_{k}^{T}\mathbf{D}_{N}^{-1}(w)\mathbf{I}_{N}\mathbf{D}_{N}(w)\boldsymbol{\phi}_{k}^{\ast}$,
independent of the candidate CFO $w$. When $w=w_{0}$, the result of
$\mathbf{D}_{N}(w_{0}-w)$ is an identity matrix, and the cost function is at a
minimum, i.e., $J(w_{0})=\sigma^{2}(N-K-L)$.

###  III-B CFO Identifiability

The CFO identifiability issue arises for the CFO estimator 9, if the following
necessary condition is fulfilled for some $w\neq w_{0}$

$\left\\{\boldsymbol{\phi}_{\widetilde{k}}^{T}\mathbf{D}_{N}(w_{0}-w)\right\\}_{\widetilde{k}\in\mathcal{N}_{\phi}}\subset\\{\boldsymbol{\phi}_{k}^{T}\\}_{k\in\mathcal{N}_{\phi}}.$
(12)

The condition in (12) for some $w\neq w_{0}$ results in multiple minimums.
Based on the definition of $\mathbf{\Phi}$ in 1, it is clear that
$\boldsymbol{\phi}_{k}$ is a chirp sequence and $\boldsymbol{\phi}_{k+\Delta
k}$ is a circularly shifted version of $\boldsymbol{\phi}_{k}$, where
$\boldsymbol{\phi}_{k+\Delta k}$ is obtained by circularly shifted
$\boldsymbol{\phi}_{k}$ down (up) by $\Delta k$ samples, with $\Delta k$ being
a non-negative (non-positive) integer. Hence, supposing that $k=K+L+\Delta k$,
the LNS of the results of $\bar{\mathbf{R}}_{yy}$ can be rewritten as

$\\{\boldsymbol{\phi}_{k}^{T}\\}_{k\in\mathcal{N}_{\phi}}=\\{\boldsymbol{\phi}_{K+L+\Delta
k}^{T}\\}_{\Delta k=1}^{N-K-L}.$ (13)

Assume a block size of $N$ with even value as an illustrative example. The
$m$th entry of $\boldsymbol{\phi}_{K+L+\Delta k}$ is given by:

$\displaystyle\left[\boldsymbol{\phi}_{K+L+\Delta k}\right]_{m}$ (14)
$\displaystyle=$
$\displaystyle\frac{e^{-\frac{j\pi}{4}}}{\sqrt{N}}e^{j\frac{\pi}{N}(K+L+\Delta
k-m)^{2}},$ $\displaystyle=$
$\displaystyle\frac{e^{-\frac{j\pi}{4}}}{\sqrt{N}}e^{j\frac{\pi}{N}(K+L-m)^{2}+j\frac{\pi}{N}(\Delta
k^{2}+2\Delta k(K+L)-2m\Delta k)},$ $\displaystyle m=1,\dots,N,\Delta
k=1,\dots,N-K-L.$

On the other hand, the sequence
$\boldsymbol{\phi}_{\widetilde{k}}^{T}\mathbf{D}_{N}(w_{0}-w)$ is also a
shifted $\boldsymbol{\phi}_{\widetilde{k}}^{T}$. Let us define
$\widetilde{k}=K+L+\Delta\widetilde{k}$, where
$\Delta\widetilde{k}=1,\dots,N-K-L$. Therefore, the $m$th entry of
$\boldsymbol{\phi}_{\widetilde{k}}^{T}\mathbf{D}_{N}(w_{0}-w)=\boldsymbol{\phi}_{K+L+\Delta\widetilde{k}}^{T}\mathbf{D}_{N}(w_{0}-w)$
is given by

$\displaystyle\left[\boldsymbol{\phi}_{K+L+\Delta\widetilde{k}}^{T}\mathbf{D}_{N}(w_{0}-w)\right]_{m}$
(15) $\displaystyle=$
$\displaystyle\frac{e^{-\frac{j\pi}{4}}}{\sqrt{N}}e^{j(\frac{\pi}{N}(K+L+\Delta\widetilde{k}-m)^{2}+(w_{0}-w)m)},$
$\displaystyle=$
$\displaystyle\frac{e^{-\frac{j\pi}{4}}}{\sqrt{N}}e^{j\frac{\pi}{N}(K+L-m)^{2}}$
$\displaystyle\times
e^{j\frac{\pi}{N}(\Delta\widetilde{k}^{2}+2\Delta\widetilde{k}(K+L)-2m\Delta\widetilde{k}+\frac{N}{\pi}(w_{0}-w)m)},$
$\displaystyle m=1,\dots,N,\Delta\widetilde{k}=1,\dots,N-K-L.$

Comparing 14 and 15, we observe that
$\boldsymbol{\phi}_{K+L+\Delta\widetilde{k}}^{T}\mathbf{D}_{N}(w_{0}-w)$ would
belong to the LNS of $\bar{\mathbf{R}}_{yy}$ if and only if the two following
conditions are fulfilled at the same time

$\displaystyle\frac{\pi}{N}(\Delta\widetilde{k}^{2}+2\Delta\widetilde{k}(K+L))\,\,{\rm
mod}\,\,2\pi$ $\displaystyle=$ $\displaystyle\frac{\pi}{N}(\Delta
k^{2}+2\Delta k(K+L))\,\,{\rm mod}\,\,2\pi,$ (16)

and

$\displaystyle\frac{2\pi}{N}m\Delta k\,\,{\rm mod}\,\,2\pi$ $\displaystyle=$
$\displaystyle(\frac{2\pi}{N}m\Delta\widetilde{k}+(w_{0}-w)m)\,\,{\rm
mod}\,\,2\pi$ (17)

for $\forall\Delta k$ and $\forall\Delta\widetilde{k}\in\\{1,\dots,N-K-L\\}$.
Recall that $(w_{0}-w)\in(-2\pi,2\pi)$. It is easy to verify that the
aforementioned two conditions can be satisfied simultaneously if and only if
$\Delta k=\Delta\widetilde{k}$, and $(w_{0}-w)\,\,{\rm mod}\,\,2\pi=0$, i.e.,
$w_{0}=w$. Therefore, the CFO estimator employing the cost function 10 has a
unique minimum and the CFO identifiability issue is resolved.

We conclude with the following proposition.

Proposition 1: For the cost function $J(w)$ in 10 to have a unique minimum,
the insertion of at least $L+1$ consecutive null subchirps guarantees the CFO
acqusition range $[-\pi,\pi)$ for an OCDM system under a multipath channel of
order up to $L$.

###  III-C Error Performance

Next, we explore the performance of the proposed OCDM transmission scheme with
consecutive null subchirps (named OCDM-NSC) w.r.t. multipath diversity. A well
known fact is that null subcarriers do not affect BER performance in OFDM
w.r.t. multipath diversity. This is to say that the multipath diversity is
always $1$ for OFDM systems due to the inherent diagonal structure of the
equivalent channel. On the other hand, null subchirps affect the OCDM
performance w.r.t. multipath diversity differently. The following proposition
is herein.

Proposition 2: Suppose that the number of consecutive null subchirps is
greater than or equal to the channel order $L$, i.e., $N-K\geq L$, the OCDM-
NSC scheme achieves full multipath diversity by LEs and MLE.

Proof: There is a link between an OCDM with null subchirps and a ZP-only
single carrier (SC) transmission scheme. It turns out that the demodulatd OCDM
block with null subchirps also assumes the same mathematical form as the ZP-
only SC block under multipath channels:

$\displaystyle\mathbf{z}(i)=$ $\displaystyle\mathbf{\Phi}\mathbf{y}(i)$
$\displaystyle=$
$\displaystyle\mathbf{\Phi}\mathbf{H}\mathbf{\Phi}^{\mathcal{H}}\mathbf{T}_{\rm
zp}\mathbf{s}(i)+\mathbf{n}(i)$ $\displaystyle=$
$\displaystyle\mathbf{H}\mathbf{T}_{\rm zp}\mathbf{s}(i)+\mathbf{n}(i),$

where $\mathbf{H}=\mathbf{\Phi}\mathbf{H}\mathbf{\Phi}^{\mathcal{H}}$ by the
property of the Fresnel matrix and the circulant matrix. It is observed in
[14, 15, 16] that ZP-only SC transmission enables full multipath diversity,
which can be collected by LEs and MLE. Therefore, the conclusion applies to
OCDM-NSC.

##  IV Simulations

In this section, we illustrate the performance results of the proposed method
through simulations. In particular, the covariance matrix $\mathbf{R}_{yy}$ is
calculated empirically across $N_{b}$ blocks as:

$\hat{\mathbf{R}}_{yy}=\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}\mathbf{y}(i)\mathbf{y}^{\mathcal{H}}(i).$
(18)

For all the simulations, we set that $N_{b}=1000$. The SNR is defined as
$E_{s}/\sigma^{2}$.

To evaluate the performance of the CFO estimation and verify its
identifiability, we define the mean square error (MSE) as
$\sum_{m=1}^{M}(w_{0}^{(m)}-\hat{w}_{0}^{(m)})^{2}/M$, where $M$ is the number
of Monte Carlo (MC) runs, and $w_{0}^{(m)}$ and $\hat{w}_{0}^{(m)}$ are the
true and estimated CFO for the $m$th MC run, respectively. Multipath channels
of order $L=2$ with Rayleigh fading coefficients are employed. The block size
is $N=16$. We set $K=12$ such that $N-K=4>L+1$ null subchirps for the proposed
CFO estimator. As our proposed CFO estimator is pilot-free, we adopt the CP-
based CFO estimator [5] for comparison. The performance of the CP-based CFO
estimator depends on the length of the excessive CP and the number of
transmission blocks. Here we use an excessive CP of length $4>(L+1)$ and
$1000$ transmission blocks for the CP-based CFO estimator for a fair
comparison. Although the CP-based CFO estimator [5] is in a closed form, it
faces the identifiability issue when $|w_{0}|>\pi/N$. The comparison results
are shown in Fig. 1, where three cases of $w_{0}\in[-0.05\pi,0.05\pi)$,
$w_{0}\in[-0.1\pi,0.1\pi)$, and $w_{0}\in[-\pi,\pi)$ are presented,
respectively. As one would expect, the CP-based CFO estimator outperforms the
proposed one, when $-\pi/N\leq w_{0}<\pi/N$ in the first case. In contrast,
when $|w_{0}|>\pi/N$, the CFO identifiability issue arises for the CP-based
CFO estimator, and the ambiguity error dominates its performance and results
in an error floor. However, the identifiability issue is well addressed for
the proposed CFO estimator. Its performance is robust for the full range of
the CFO, i.e., $w_{0}\in[-\pi,\pi)$. It is also worth pointing out that our
proposed CFO estimator can work well in complement with the CP-based CFO
estimator. The proposed method addresses the identifiability issue and achieve
a coarse CFO compensation, and subsequently the CP-based method [5] carries
out a finer tuning. This resulting two-step CFO compensation scheme reaps both
high accuracy and robustness against a full range CFO.

Figure 1: MSEs of CFO estimators versus SNR.

Next, we simulate the BER performance of the OCDM-NSC system with the proposed
CFO estimator over multipath channels with Rayleigh fading coefficients. The
QPSK modulation is adopted. The LEs specified in 3 and 4 are applied,
respectively. The true CFO $w_{0}$ is randomly generated following a uniform
distribution in the range of $[-\pi,\pi)$ for each MC run. The simulation
parameters are set as $N=16$, $K=12$, and $L=2$. For a fair comparison, the
OFDM performance using the same number of null subcarriers and with distinct
spacing [10], and the performance of the CP-based CFO estimation method [5]
are also plotted in Fig. 2. Evidently from Fig. 2, all three equalizers for
OCDM-NSC collect full diversity. This is a direct result of the tall Toeplitz
structure of the equivalent channel for OCDM-NSC, i.e.,
$\mathbf{H}\mathbf{T}_{\rm zp}$. It provides a better channel matrix
condition, and results a better BER performance [15, 16]. With this setup,
OFDM asymptotically have unit diversity, whereas the performance of the CP-
based CFO estimation method suffers because of the identifiability ambiguity.
In this regard, we have testified that the proposed CFO estimator guarantees
the CFO identifiability. Meanwhile, the null subchirps enable the multipath
diversity of the OCDM-NSC system, and MLE and LEs collect the diversity gain.

Figure 2: BER performance with the proposed CFO estimator.

##  V Conclusion

This paper investigates and solves the CFO identifiability problem in OCDM
systems. We propose to insert consecutive null subchirps to facilitate the CFO
estimation. These null subchirps restore CFO identifiability by creating
channel independent null subspace to overcome adverse channel conditions. A
CFO estimator is accordingly proposed to achieve a full acquisition range. It
has also been demonstrated that the OCDM system with consecutive null
subchirps enables better performance than its plain counterpart. Finally, CFO
identifiability for the proposed estimator is validated through simulations.
Performance comparisons are made against other methods to show the advantages
of the proposed OCDM-NSC.

## References

  * [1] X. Ouyang and J. Zhao, “Orthogonal chirp division multiplexing,” _IEEE Transactions on Communications_ , vol. 64, no. 9, pp. 3946–3957, 2016.
  * [2] X. Ouyang, O. A. Dobre, Y. L. Guan, and J. Zhao, “Chirp spread spectrum toward the nyquist signaling rate—orthogonality condition and applications,” _ieee signal processing letters_ , vol. 24, no. 10, pp. 1488–1492, 2017.
  * [3] M. S. Omar and X. Ma, “Designing ocdm-based multi-user transmissions,” in _2019 IEEE Global Communications Conference (GLOBECOM)_. IEEE, 2019, pp. 1–6.
  * [4] X. Ouyang, C. Antony, G. Talli, and P. D. Townsend, “Robust channel estimation for coherent optical orthogonal chirp-division multiplexing with pulse compression and noise rejection,” _Journal of Lightwave Technology_ , vol. 36, no. 23, pp. 5600–5610, 2018.
  * [5] R. Zhang, Y. Wang, and X. Ma, “Channel estimation for ocdm transmissions with carrier frequency offset,” _IEEE Wireless Communications Letters_ , vol. 11, no. 3, pp. 483–487, 2021.
  * [6] M. L. de Filomeno, L. G. de Oliveira, Â. Camponogara, A. Diewald, T. Zwick, M. L. de Campos, and M. V. Ribeiro, “Joint channel estimation and schmidl & cox synchronization for ocdm-based systems,” _IEEE Communications Letters_ , vol. 26, no. 8, pp. 1878–1882, 2022.
  * [7] L. G. de Oliveira, B. Nuss, M. B. Alabd, A. Diewald, Y. Li, L. Gehre, X. Long, T. Antes, J. Galinsky, and T. Zwick, “Discrete-fresnel domain channel estimation in ocdm-based radar systems,” _IEEE Transactions on Microwave Theory and Techniques_ , 2022.
  * [8] M. S. Omar and X. Ma, “Performance analysis of ocdm for wireless communications,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 7, pp. 4032–4043, 2021.
  * [9] V. K. Trivedi and P. Kumar, “Low complexity interference compensation for dfrft-based ofdm system with cfo,” _IET Communications_ , vol. 14, no. 14, pp. 2270–2281, 2020.
  * [10] X. Ma, C. Tepedelenlioglu, G. B. Giannakis, and S. Barbarossa, “Non-data-aided carrier offset estimators for ofdm with null subcarriers: identifiability, algorithms, and performance,” _IEEE Journal on selected areas in communications_ , vol. 19, no. 12, pp. 2504–2515, 2001.
  * [11] M. Ghogho, A. Swami, and G. B. Giannakis, “Optimized null-subcarrier selection for cfo estimation in ofdm over frequency-selective fading channels,” in _GLOBECOM’01. IEEE Global Telecommunications Conference (Cat. No. 01CH37270)_ , vol. 1. IEEE, 2001, pp. 202–206.
  * [12] F. Gao and A. Nallanathan, “Identifiability of data-aided carrier-frequency offset estimation over frequency selective channels,” _IEEE transactions on signal processing_ , vol. 54, no. 9, pp. 3653–3657, 2006.
  * [13] F. Gao, T. Cui, and A. Nallanathan, “Scattered pilots and virtual carriers based frequency offset tracking for ofdm systems: algorithms, identifiability, and performance analysis,” _IEEE Transactions on Communications_ , vol. 56, no. 4, pp. 619–629, 2008.
  * [14] Z. Wang and G. B. Giannakis, “Wireless multicarrier communications,” _IEEE signal processing magazine_ , vol. 17, no. 3, pp. 29–48, 2000.
  * [15] C. Tepedelenlioglu and Q. Ma, “On the performance of linear equalizers for block transmission systems,” in _GLOBECOM’05. IEEE Global Telecommunications Conference, 2005._ , vol. 6. IEEE, 2005, pp. 5–pp.
  * [16] X. Ma and W. Zhang, “Fundamental limits of linear equalizers: diversity, capacity, and complexity,” _IEEE Transactions on Information Theory_ , vol. 54, no. 8, pp. 3442–3456, 2008.