Jakob Assländer

*Jakob Assländer, Center for Biomedical Imaging, Department of Radiology, New York University Grossman School of Medicine, 650 1st Avenue, Room 218, New York, NY 10016, USA.

grant NIH/NIBIB R21 EB027241 and was performed under the rubric of the Center
for Advanced Imaging Innovation and Research, a NIBIB Biomedical Technology
Resource Center (NIH P41 EB017183).

# Generalized Bloch model: a theory for pulsed magnetization transfer

Jakob Assländer  Cem Gultekin  Sebastian Flassbeck  Steffen J Glaser  Daniel K
Sodickson  Center for Biomedical Imaging, Dept. of Radiology, New York
University Grossman School of Medicine, NY, USA  Center for Advanced Imaging
Innovation and Research (CAI2R), Dept. of Radiology, New York University
Grossman School of Medicine, NY, USA  Courant Institute of Mathematical
Sciences, New York University, NY, USA  Department of Chemistry, Technische
Universität München, Germany<EMAIL_ADDRESS>
###### Abstract

[Abstract] Purpose: The paper introduces a classical model to describe the
dynamics of large spin-1/2 ensembles associated with nuclei bound in large
molecule structures, commonly referred to as the semi-solid spin pool, and
their magnetization transfer (MT) to spins of nuclei in water.

Theory and Methods: Like quantum-mechanical descriptions of spin dynamics and
like the original Bloch equations, but unlike existing MT models, the proposed
model is based on the algebra of angular momentum in the sense that it
explicitly models the rotations induced by radio-frequency (RF) pulses. It
generalizes the original Bloch model to non-exponential decays, which are,
e.g., observed for semi-solid spin pools. The combination of rotations with
non-exponential decays is facilitated by describing the latter as Green’s
functions, comprised in an integro-differential equation.

Results: Our model describes the data of an inversion-recovery magnetization-
transfer experiment with varying durations of the inversion pulse
substantially better than established models. We made this observation for all
measured data, but in particular for pulse durations small than 300$\upmu$s.
Furthermore, we provide a linear approximation of the generalized Bloch model
that reduces the simulation time by approximately a factor 15,000, enabling
simulation of the spin dynamics caused by a rectangular RF-pulse in roughly
2$\upmu$s.

Conclusion: The proposed theory unifies the original Bloch model, Henkelman’s
steady-state theory for magnetization transfer, and the commonly assumed
rotation induced by hard pulses (i.e., strong and infinitesimally short
applications of RF fields) and describes experimental data better than
previous models.

###### keywords:

quantitative magnetization transfer, MT, qMT, chemical exchange, quantitative
MRI, parameter mapping

††articletype: Pre-Print - submitted to Magnetic Resonance in
Medicine††footnotetext: Word Count: 5169

## Introduction

Commonly, magnetization transfer 1 is described by a 2-pool model 2 that
distinguishes between protons bound in water—the so-called free pool—and
protons bound in macromolecules, such as proteins or lipids—the so-called
semi-solid pool. The transversal magnetization of the semi-solid pool decays
rather quickly ($T_{2}^{s}\approx 10\upmu$s). Therefore, we commonly do not
observe it directly. Nevertheless, the signal or contrast in virtually all MR
images is affected by these molecules 3, 4 due to the indirect effect of
magnetization transfer,1 i.e. the transfer of z-magnetization between the two
pools.

The original Henkelman model 2 describes the spin-dynamics of these two pools
jointly with the Bloch-McConnell equations.5 During a continuous radio
frequency wave, a steady state evolves, which allows for solving the system of
differential equations. One of the key contributions of Henkelman et al. 2 was
to reformulate the steady-state solution of the Bloch-McConnell equations to
isolate two Lorentzian lineshapes that describe the saturation of respective
spin pools by the continuous RF wave. This mathematical trick views the
experiment from the perspective of absorption NMR and, in line with the theory
of absorption NMR, Henkelman et al. replaced the Lorentzian with a Gaussian
lineshape to more accurately describe the signal of an agar phantom,2 and with
a super-Lorentzian lineshape to describe the signal of brain white matter.6

In gases and liquids, fast molecular motion allows one to model nuclear
dipole-dipole interactions as random magnetic field fluctuations 7 and to
capture their macroscopic effect by two relaxation rates, $R_{1}=1/T_{1}$ and
$R_{2}=1/T_{2}$. This theory forms the basis for the original Bloch equations
8 and their variants, such as the Bloch-McConnell equations.5 Equivalently,
the fast molecular motion gives rise to the Lorentzian lineshape in
Henkelman’s RF-absorption picture.

The molecular motion of large proteins or lipids is, however, substantially
slower and, as a consequence, dipole-dipole interactions cannot be
approximated by random fields. Thus, the macroscopic magnetization dynamics
cannot be captured by the traditional Bloch equations or a corresponding
Lorentzian lineshape. In cases where the Bloch equations fail, one can usually
rely on quantum mechanics to provide an accurate model. However, the large
number of interacting nuclei in each molecule and the large number of
different macromolecules present in biological tissue render a full quantum
mechanical treatment infeasible for the MT effect.

Henkelman et al. were able to avoid a quantum mechanical description by
focusing on a steady-state condition in the presence of continuous RF waves—a
condition in which absorption NMR theory applies. Clinical MRI, however, is
performed almost exclusively with RF-pulses9 instead of continuous RF waves,
due to safety limitations associated with RF power deposition, and due to
practical time constraints. In recognition of the need for a classical pulsed-
MT model, Graham et al.10 and Sled and Pike11 modified Henkelman’s theory and
demonstrated good agreement with experimental observations for off-resonant
saturation pulses with a duration of several milliseconds.

For short RF-pulses, however, these models fail to describe the spin dynamics
accurately. This becomes apparent when approaching the extreme case of a hard
pulse, i.e. a strong and infinitesimally short application of RF field.
Manning et al.12 demonstrated experimentally that a hard inversion pulse does
indeed invert the semi-solid pool, which stands in contrast to existing MT
models that predict a saturation of the semi-solid pool to zero in this limit.

Here, we propose a model that improves existing MT models by the following
properties:

  * •

It models the limit of a hard pulse correctly.

  * •

It interpolates between the extremes of a hard pulse and a continuous wave,
and the modeled spin dynamics is in line first experimental validations.

  * •

It is equivalent to the Bloch equations when assuming a Lorentzian lineshape.

It is a classical theory that generalizes the Bloch equations to arbitrary
lineshapes to capture the macroscopic effect of dissipative dipole-dipole
interactions in large molecule structures.

## Theory

### Isolated semi-solid pool

We first derive the generalized Bloch equations for an isolated semi-solid
pool before incorporating magnetization transfer to and from the free pool.

Expressing the Bloch equations8 as

$\displaystyle(\partial_{t}+R_{2}^{s})x^{s}(t)$
$\displaystyle=\omega_{y}(t)z^{s}(t)$ (1)
$\displaystyle(\partial_{t}+R_{2}^{s})y^{s}(t)$
$\displaystyle=-\omega_{x}(t)z^{s}(t)$ (2)
$\displaystyle(\partial_{t}+R_{1}^{s})z^{s}(t)$
$\displaystyle=-\omega_{y}(t)x^{s}(t)+\omega_{x}(t)y^{s}(t)+R_{1}^{s}$ (3)

isolates the linear differential operators $\partial_{t}+R_{1,2}^{s}$, which
are composed of the partial derivative wrt. time and relaxation. Here,
$\omega_{x,y}$ are the Rabi frequencies 13 in both dimensions. Without loss of
generality, we have described the spin system in a frame of reference that
rotates at the Larmor frequency, hence $\omega_{z}=0$. The three spatial
components of the magnetization are denoted by $x^{s}(t),y^{s}(t),z^{s}(t)$
and are normalized by setting the thermal equilibrium magnetization to
$z^{s}(0)=1$.

Eqs. (1) and (2) are formally solved as follows:

$\displaystyle x^{s}(t)$
$\displaystyle=\int_{0}^{t}G(t,\tau)\omega_{y}(\tau)z^{s}(\tau)d\tau$ (4)
$\displaystyle y^{s}(t)$
$\displaystyle=-\int_{0}^{t}G(t,\tau)\omega_{x}(\tau)z^{s}(\tau)d\tau,$ (5)

where

$G(t,\tau)=\exp(-R_{2}^{s}(t-\tau))\;\;\forall\;\;t\geq\tau$ (6)

is the Green’s function of the operator $\partial_{t}+R_{2}^{s}$ when assuming
thermal equilibrium as initial conditions ($x^{s}(0)=y^{s}(0)=0$,
$z^{s}(0)=1$). Physically, we can interpret this solution by splitting the
radio-frequency field $\omega_{x,y}$ into infinitesimally short hard pulses:
Each hard pulse at the time $\tau$ excites spin coherence or transversal
magnetization, which, thereafter, decays freely for the time $t-\tau$ before
it is observed at the time $t$. Since $G(t,\tau)=G(t-\tau)$ is only a function
of the time between excitation and observation, Eqs. (4) and (5) describe a
convolution of the Green’s function (Fig. 1)—which models an FID—with the
time-varying RF-field and the evolution of the $z^{s}$-magnetization.

Figure 1: Green’s Functions. They are given by the Fourier transform of
respective lineshape and are equivalent to the free induction decay of the
transversal magnetization after excitation with a hard, i.e. infinitesimally
short, RF-pulse (disregarding inhomogeneous broadening, i.e. assuming
$T_{2}^{*,s}=T_{2}^{s}$). Despite having the same $T_{2}^{s}=10\upmu$s, the
decay characteristics differ substantially: while the magnetization decays
quickly when assuming a Lorentzian or Gaussian lineshape, more than 1% of the
magnetization remains after 1ms when assuming a super-Lorentzian lineshape.

For hard pulses with an infinitesimally small flip angle, Ernst et al. showed
that the FID is equivalent to the Fourier transform of the spectral absorption
lineshape.14, 15 Therefore, we make the ansatz that we can mimic Henkelman’s
approach and replace the Green’s function of Eq. (6) with the Fourier
transform of other lineshapes, such as the Gaussian curve

$G(t,\tau)=\exp(-{R_{2}^{s}}^{2}(t-\tau)^{2}/2)),$ (7)

which describes the FID of gels well.2 Inserting this Green’s function into
Eqs. (4) and (5) facilitates the description of the spin dynamics in gels
during a finite RF-pulse.

As shown by Morrison et al.,6 the absorption lineshape of brain white matter
is well-described by a super-Lorentzian lineshape.16 Hence, its Fourier
transform provides us with the Green’s function

$G(t,\tau)=\int_{0}^{1}\exp\left(-{R_{2}^{s}}^{2}(t-\tau)^{2}\cdot\frac{(3\zeta^{2}-1)^{2}}{8}\right)d\zeta,$
(8)

which approaches the function
$G(t,\tau)=\sqrt{\frac{2\pi}{3R_{2}^{s}(t-\tau)}}$ for $t-\tau\gg T_{2}^{s}$.
All three Green’s functions are depicted in Fig. 1.

The longitudinal magnetization’s return to thermal equilibrium is well-
described by a mono-exponential function.17 Hence, Eqs. (3)-(5), which we dub
the generalized Bloch model, describe the dynamics of an isolated semi-solid
pool. As the transversal magnetization decays before we can measure it with
common clinical pulse sequences, one is often interested only in its implicit
effect on the longitudinal magnetization. Hence, we can insert Eqs. (4) and
(5) into Eq. (3) to yield a single equation that describes the dynamics of the
longitudinal magnetization:

$\begin{split}\partial_{t}z^{s}(t)=-\omega_{y}(t)\int_{0}^{t}G(t,\tau)\omega_{y}(\tau)z^{s}(\tau)d\tau\\\
-\omega_{x}(t)\int_{0}^{t}G(t,\tau)\omega_{x}(\tau)z^{s}(\tau)d\tau\\\
+R_{1}^{s}(1-z^{s}(t)).\end{split}$ (9)

This is an integro-differential equation, which can be solved numerically.

Even though Eq. (9) describes the dynamics of $z$-magnetization without
explicitly denoting transversal magnetization, we note that transverse
magnetization components are implicitly baked into this equation, and the
multiplication of $z^{s}(\tau)$ with $\omega_{x,y}(\tau)$ as well as with
$\omega_{x,y}(t)$ reflects the angular momentum algebra.

### Magnetization transfer to and from the free pool

So far, we have assumed an isolated semi-solid spin pool. In order to describe
magnetization transfer, we incorporate the generalized Bloch equations into a
2-pool model. We choose this model over the more accurate 4-pool model
described by Manning et al.12 because it is difficult to determine the large
number of parameters of the 4-pool model in vivo due to experimental
constraints. However, the generalized Bloch model can be incorporated in
models with an arbitrary number of pools.

We can describe the three spatial components of the free pool’s magnetization
($x^{f}$, $y^{f}$, and $z^{f}$) by the original Bloch equations due to the
fast molecular motion associated with the free pool. Since the transversal
magnetization of the semi-solid pool is short-lived ($R_{2}^{s}\gg R_{x}$,
where $R_{x}$ denotes the exchange rate between the two pools), we can neglect
coherence transfer between the two pools and it is sufficient to account for
the dynamics of the z-component $z^{s}$.2 Given a thermal equilibrium
magnetization of the free and semi-solid pool of $m_{0}^{f}$ and $m_{0}^{s}$,
respectively, and normalizing their sum as $m_{0}^{f}+m_{0}^{s}=1$, we can
describe the dynamics of the coupled spin system by the generalized Bloch-
McConnell equations:

$\displaystyle(\partial_{t}+R_{2}^{f})x^{f}(t)=$
$\displaystyle\;\omega_{y}(t)z^{f}(t)$ (10)
$\displaystyle(\partial_{t}+R_{2}^{f})y^{f}(t)=$
$\displaystyle-\omega_{x}(t)z^{f}(t)$ (11)
$\displaystyle(\partial_{t}+R_{1}^{f})z^{f}(t)=$
$\displaystyle-\omega_{y}(t)x^{f}(t)+\omega_{x}(t)y^{f}(t)$ (12)
$\displaystyle+m_{0}^{f}R_{1}^{f}+R_{x}(m_{0}^{f}z^{s}(t)-m_{0}^{s}z^{f}(t))$
$\displaystyle(\partial_{t}+R_{1}^{s})z^{s}(t)=$
$\displaystyle-\omega_{y}(t)\int_{0}^{t}G(t,\tau)\omega_{y}(\tau)z^{s}(\tau)d\tau$
(13)
$\displaystyle-\omega_{x}(t)\int_{0}^{t}G(t,\tau)\omega_{x}(\tau)z^{s}(\tau)d\tau$
$\displaystyle+m_{0}^{s}R_{1}^{s}+R_{x}(m_{0}^{s}z^{f}(t)-m_{0}^{f}z^{s}(t)).$

## Methods

### Simulations

The simulations presented here were designed to highlight the generalized
Bloch model’s ability to interpolate between the well-understood extrema of a
hard (i.e. infinitesimally short) RF-pulse and a continuous (i.e. infinitely
long) RF-wave, and to show agreements with and deviations from existing
models. For this purpose, we focus here on rectangular RF-pulses but note that
the generalized Bloch model can be applied to arbitrary pulse shapes. The
validity of the generalized Bloch equations for different pulse shapes is
briefly touched upon in the Discussion section.

We simulated the dynamics of an isolated semi-solid spin-pool for two
scenarios:

  * •

The spin-pool’s approach to a steady state assuming a continuous wave with an
off-resonance frequency of $\Delta/(2\pi)=100$Hz and an amplitude of
$\omega_{1}/(2\pi)=100$Hz, as well as with an amplitude of
$\omega_{1}/(2\pi)=1000$Hz.

  * •

The spin-pool’s state at the end of on-resonant RF-pulses ($\Delta=0$) with
flip angles of $\alpha=\pi/4$, $\alpha=\pi/2$, and $\alpha=\pi$ and varying
pulse durations $T_{\text{RF}}\in[0.1\upmu\text{s},0.1\text{s}]$. As we are
simulating rectangular pulses, their amplitude is given by
$\omega_{1}=\alpha/T_{\text{RF}}$.

Amplitude and frequency can be translated to the rotating frame of reference
used in Eq. (9) via the expressions $\omega_{x}=\omega_{1}\cos(\Delta t)$ and
$\omega_{y}=\omega_{1}\sin(\Delta t)$. For both of our selected scenarios we
used a subset of the models described below assuming Lorentzian, Gaussian, and
super-Lorentzian lineshapes, each with $T_{2}^{s}=10\upmu$s
($R_{2}^{s}=10^{5}/$s) and $R_{1}^{s}=1/$s.

#### Generalized Bloch model

We solved Eq. (9) with the DifferentialEquations.jl18 package. It provides a
solver for delay differential equations and we represent $z^{s}(\tau)$ as an
interpolated history function over which we numerically integrate.

#### Bloch model

In order to verify the equivalence of the generalized Bloch model using a
Lorentzian lineshape to the original Bloch equations (Eqs. (1)-(3)), we
performed Bloch simulations by denoting the original Bloch equations in matrix
form (cf. Eq. (1) in Ref. 19) and taking the matrix exponential.

#### Henkelman’s steady-state model

In order to demonstrate the equivalence of the generalized Bloch model with
Henkelman’s model when the magnetization is in a steady state, we simulated
Henkelman’s model for all three lineshapes mentioned earlier. As we here
assume an isolated semi-solid pool, Henkelman’s steady-state solution (Eq. (8)
in Ref. 2) simplifies to

$R_{1}^{s}\cdot(1-z_{\text{ss}}^{s})=R_{\text{RF}}^{s}z_{\text{ss}}^{s},$ (14)

where $z_{\text{ss}}^{s}$ denotes the steady-state magnetization. Eq. (14)
describes the equilibrium of $T_{1}$ relaxation and RF-saturation with the
rate

$R_{\text{RF}}^{s}=\pi\omega_{1}^{2}g(\Delta,T_{2}^{s}),$ (15)

where $g(\Delta,T_{2}^{s})$ denotes the spectral absorption lineshape (Eqs.
(A.25)-(A.27) in the appendix).

#### Graham’s spectral model

Inspired by Henkelman’s steady-state model, Graham et al.10 proposed to model
the dynamics of the semi-solid spin pool with an exponential function:

$\partial_{t}z^{s}(t)=-R_{\text{RF}}^{s}z^{s}(t)+R_{1}^{s}(1-z^{s}).$ (16)

They postulated that the average saturation rate of an RF pulse can be
described by

$\langle
R_{\text{RF}}^{s}\rangle=\pi\int_{-\infty}^{+\infty}d\Delta^{\prime}S(\Delta^{\prime},\Delta,\omega_{1}(t))g(\Delta^{\prime},T_{2}^{s}),$
(17)

where $S(\Delta^{\prime},\Delta,\omega_{1}(t))$ denotes the power spectral
density, observed at the frequency $\Delta^{\prime}$, resulting from an RF-
pulse with the frequency $\Delta$.20 For the purpose of this paper, we call
the model described by Eqs. (16) and (17) Graham’s spectral model. $\langle
R_{\text{RF}}^{s}\rangle$ for on-resonant rectangular pulses and for all three
lineshapes are derived in the appendix (Eqs. (A.30), (A.31), and (A.32)).

#### Graham’s single-frequency approximation

In Ref. 10, Graham et al. analyzed off-resonant RF-pulses with durations of
several milliseconds. They observed that a Gaussian lineshape is approximately
constant throughout the spectral response of such RF-pulses (cf. Fig. A.9).
Mathematically, this is equivalent to assuming

$S(\Delta^{\prime},\Delta,\omega_{1})\approx\delta(\Delta^{\prime}-\Delta)\int_{0}^{T_{\text{RF}}}d\tau\omega_{1}^{2}(\tau).$
(18)

Here, $\delta(x)$ denotes the Dirac-delta distribution and the subsequent
integral can be derived from Parseval’s theorem. With this assumption, the
average saturation rate simplifies to

$\langle R_{\text{RF}}^{s}\rangle=\frac{\pi
g(\Delta,T_{2}^{s})}{T_{\text{RF}}}\int_{0}^{T_{\text{RF}}}d\tau\omega_{1}^{2}(\tau).$
(19)

From this average saturation rate, Graham et al. heuristically deduced a time-
dependent saturation rate of

$R_{\text{RF}}^{s}(t)=\pi\omega_{1}^{2}(t)g(\Delta,T_{2}^{s}),$ (20)

which is essentially a time-dependent version of Henkelman’s steady-state
saturation rate (Eq. (15)). Note that the last step is strictly true only for
a constant $\omega_{1}$, which is to say a rectangular RF-pulse. For the
purpose of this paper, we call the model described by Eqs. (16) and (20)
Graham’s single-frequency approximation.

Note that Graham’s single-frequency approximation is not valid for a super-
Lorentzian lineshape and on-resonant RF-pulses due to the divergence of the
lineshape. A more detailed spectral analysis of this approximation can be
found in the appendix (Fig. A.9).

#### Sled’s model

Sled and Pike11 used the same exponential model as Graham and proposed to
describe the time-dependent saturation rate by a convolution of the squared
RF-pulse amplitude with the impulse response $G(t-\tau)$:

$\partial_{t}z^{s}(t)=\left(-\pi\int_{0}^{t}G(t-\tau)\omega_{1}(\tau)^{2}d\tau\right)z^{s}(t)+R_{1}^{s}(1-z^{s}).$
(21)

The impulse response used by Sled and Pike is equivalent to the Green’s
functions derived here, and is indicated using the same notation.

### Experiments

In order to verify the generalized Bloch equations, we performed inversion
recovery experiments21 with varying inversion-pulse durations on a 1T NMR
spectrometer (Spinsolve 43MHz, Magritek, Wellington, New Zealand).

#### NMR samples

We built an NMR sample of thermally cross-linked bovine serum albumin (BSA),
where the BSA makes up 15% of the overall sample weight. We mixed the BSA
powder with distilled water (weight ratio 15 to 85) and stirred it at 30∘C
until the BSA was fully dissolved. We filled an NMR tube with the resulting
liquid and thermally cross-linked it in a water bath at approximately 90∘C for
10 minutes.

For comparison, we also built an NMR sample consisting of distilled water
doped with 0.15mM MnCl2.

#### Pulse sequence

We implemented a selective inversion recovery (SIR) experiment similar to the
one proposed by Gochberg et. al.21 Starting from thermal equilibrium (ensured
by a repetition time of 30s), the free pool is inverted with a rectangular
$\pi$-pulse. After an inversion time $T_{\text{i}}$, we excite the
magnetization with a narrow-band composite $\pi/2$-pulse,22 consisting of 12
$\pi$-pulses with the RF-phases (104.5 + 104.5)∘, (104.5 + 255.5)∘, 104.5∘,
(104.5 + 255.5)∘, (104.5 + 104.5)∘, (255.5 + 104.5)∘, (255.5 + 255.5)∘,
255.5∘, (255.5 + 255.5)∘, (255.5 + 104.5)∘, 104.5∘, and 255.5∘, followed by a
$\pi/2$-pulse with the zero phase. The goal of this composite pulse is to
excite only magnetization with a well-calibrated $B_{1}^{+}$-field or,
equivalently, only magnetization that was, to good approximation, fully
inverted by the rectangular $\pi$-pulse.

We use phase cycling to select the FID of the excitation pulse while canceling
out the FID of the (imperfect) inversion pulse, as well as the spin echo of
their combination. The excitation profile of the composite pulse has a small
spurious imaginary part in addition to the desired large real part.22 We
remove this imaginary part with phase cycling as well. Together, this requires
8 repetitions of each experiment. We acquire the FID for 1.6s with a dwell
time of 100$\upmu$s.

We measured 200 data points overall, with 20 different inversion times
$T_{\text{i}}\in[3\text{ms},5\text{s}]$, scaled exponentially, and 10
different durations of the inversion pulse:
$T_{\text{RF}}\in\\{22.8,45.6,114,228,342,456,570,684,798,912\\}\upmu$s.

#### Data analysis

We performed the following analyses of the inversion-recovery (IR) NMR data.
All fits were performed with the Levenberg-Marquardt algorithm and we
calculated an uncertainty from the Jacobian matrix at the point of
convergence. For more details, refer to the published source code.

First, we estimated $T_{2}^{*,f}$ by fitting a mono-exponential decay curve to
the FID of the acquisition with $T_{\text{RF}}=22.8\upmu$s and
$T_{\text{i}}=5$s.

For the remaining analyses, we considered only the first data point of each
FID. Plotting these data points as a function of their inversion time
$T_{\text{i}}$ results in one inversion recovery curve for each
$T_{\text{RF}}$. The phase of the signal is overall very stable, which allowed
us to multiply each inversion recovery curve by the complex conjugate phase
factor of the longest $T_{\text{i}}$ data point and to disregard the imaginary
part, which contains mostly noise.

We analyzed the resulting IR curves for each $T_{\text{RF}}$ separately by
fitting to a mono-exponential function. Further, we performed a global fit on
all data for all $T_{\text{i}}$ and $T_{\text{RF}}$. We fitted the coupled
spin-pool model (Eqs. (10)-(13)), which describes the free pool with the
original Bloch model and the semi-solid pool with the generalized Bloch model,
assuming a super-Lorentzian lineshape. The super-Lorentzian lineshape for the
BSA sample was confirmed with a continuous-wave saturation experiment (the
analysis can be found on
https://jakobasslaender.github.io/MRIgeneralizedBloch.jl/v0.2.11/build_literate/Analyze_NMR_PreSat_Data).
We repeated the fits with the same coupled spin-pool model while replacing the
description of the semi-solid pool with Graham’s spectral model (Eqs.
(16),(17)) and Sled’s model (Eq. (21)), respectively.

All considered coupled-spin models comprise 7 parameters: an overall scaling
factor, an inversion efficiency of the free pool, defined by
$-z^{f}(0)/z^{f}(+\infty)$, the macromolecular pool fraction $m_{0}^{s}$, the
longitudinal relaxation rate, where we assume $R_{1}=R_{1}^{f}=R_{1}^{s}$ as
proposed in the literature,21 $T_{2}^{f}$, $T_{2}^{s}$, and the exchange rate
$R_{x}$. We assumed $T_{2}^{f}\approx T_{2}^{*,f}$ and extracted the latter
from the FID for the following reason: a finite $T_{2}^{f}$-value causes a
reduction of the inversion efficiency with increasing $T_{\text{RF}}$ and, in
the presence of off-resonance, $T_{2}^{*,f}$ models this effect better than
$T_{2}^{f}$. We note that the effect is small regardless, as
$\\{T_{2}^{f},T_{2}^{*,f}\\}\gg T_{\text{RF}}$.

To assess how well the models describe the observed spin dynamics, we
calculate residuals: For each model and each $T_{\text{RF}}$ separately, we
calculated the $\ell_{2}$-norm of the difference between measured and modeled
signal and normalized it by the $\ell_{2}$-norm of the measured signal. First,
we calculate the residuals of the actual fits, i.e. by simulating the signal
with each model and the biophysical parameters that were estimated with
respective model. For Graham’s and Sled’s models’, this has, however, the
disadvantage that the poor description at short $T_{\text{RF}}$ negatively
affects the residuals at long $T_{\text{RF}}$. To overcomes this problem, we
additionally calculate the residuals between the measured signal and the one
simulated with each model, but with the biophysical parameters that were
estimated by fitting the generalized Bloch model.

### Linear approximation of the generalized Bloch equations

Solving Eqs. (10)-(13) is numerically challenging. For scenarios that require
the simulation of many RF-pulses, as well as a repeated simulation of these
RF-pulses (e.g., in a fitting routine for many voxels), the simulation time
becomes impractically long. This problem can be overcome with a linear
approximation of the generalized Bloch model. We hypothesize that we can
approximate the spin dynamics using a simplified Bloch-McConnell approach 5

$\partial_{t}\begin{pmatrix}x^{f}\\\ y^{f}\\\ z^{f}\\\ x^{s}\\\ z^{s}\\\
1\end{pmatrix}=\begin{pmatrix}-R_{2}^{f}&-\omega_{z}&\omega_{y}&0&0&0\\\
\omega_{z}&-R_{2}^{f}&0&0&0&0\\\
-\omega_{y}&0&-R_{1}^{f}-R_{\text{x}}m_{0}^{s}&0&R_{\text{x}}m_{0}^{f}&m_{0}^{f}R_{1}^{f}\\\
0&0&0&-R_{2}^{s,l}(T_{\text{RF}}/T_{2}^{s},\alpha)&\omega_{y}&0\\\
0&0&R_{\text{x}}m_{0}^{s}&-\omega_{y}&-R_{1}^{s}-R_{\text{x}}m_{0}^{f}&m_{0}^{s}R_{1}^{s}\\\
0&0&0&0&0&0\end{pmatrix}\begin{pmatrix}x^{f}\\\ y^{f}\\\ z^{f}\\\ x^{s}\\\
z^{s}\\\ 1\end{pmatrix}.$ (22)

As described above, we account only for exchange of the $z$-components as
$R_{2}^{s}\gg R_{\text{x}}$. Furthermore, we disregard the $y^{s}$-component
as $R_{2}^{s}\gg\omega_{z}$ in the rotating frame, and we use RF-pulses
aligned by convention along the $y$-axis.

Given this Bloch-McConnell model, the task at hand is to determine a
linearized $R_{2}^{s,l}$ that minimizes the approximation error. For this
purpose, we can solve the full generalized Bloch-McConnell model (Eqs.
(10)-(13)) and numerically search for the $R_{2}^{s,l}$ that best approximates
the magnetization vector at the end of the pulse. In general, this
magnetization vector depends on the eight parameters $\omega_{y}$,
$\omega_{z}$, $m_{0}^{s}$, $R_{1}^{f}$, $R_{1}^{s}$, $R_{2}^{f}$, $R_{2}^{s}$,
and $R_{\text{x}}$, as well as the initial magnetization, rendering a pre-
computation of $R_{2}^{s,l}$ over this entire space infeasible. Fortunately,
the substantial difference in time scales
($R_{2}^{s}\gg\\{R_{\text{x}},R_{1}^{f},R_{1}^{s}\\}$) decouples the different
processes. For the approximation of $R_{2}^{s,l}$, this allows us to treat the
semi-solid pool as isolated (Eq. (9)) and to assume $R_{1}^{s}=0$. Under these
assumptions, the magnetization at the end of an RF-pulse depends only on
$T_{2}^{s}$, $\omega_{y}$, and $T_{\text{RF}}$. In fact, it only depends on
the ratio $T_{\text{RF}}/T_{2}^{s}$ and the product $\omega_{y}T_{2}^{s}$ or,
equivalently, on $T_{\text{RF}}/T_{2}^{s}$ and the flip angle $\alpha$. This
simplification results from a transformation of Eq. (9) to a dimensionless
space, i.e. it results from substituting $\tilde{t}=t/T_{\text{RF}}$ and
$\tilde{\tau}=\tau/T_{\text{RF}}$ and integrating the differential equation
from $0$ to $T_{\text{RF}}$ while assuming a constant $\omega_{y}$. In
summary, these approximations and reformulations reduce the problem to a two-
dimensional space and we can compare the $z^{s}$ at the end of the RF-pulse,
as modeled by Eq. (9), to the one modeled by

$\partial_{t}\begin{pmatrix}x^{s}\\\
z^{s}\end{pmatrix}=\begin{pmatrix}-R_{2}^{s,l}(\frac{T_{\text{RF}}}{T_{2}^{s}},\alpha)&\omega_{y}\\\
-\omega_{y}&0\\\ \end{pmatrix}\begin{pmatrix}x^{s}\\\ z^{s}\end{pmatrix}.$
(23)

With $T_{\text{R}}\gg T_{2}^{s}$, $x^{s}$ vanishes between RF-pulses. Hence,
we can assume $x^{s}(0)=0$ at the beginning of each RF-pulse and we can,
further, neglect $x^{s}$ at the end of the pulse.

These simplifications render a pre-computation of $R_{2}^{s,l}$ feasible. We
solved Eq. (9) for 4096 combinations of $T_{\text{RF}}/T_{2}^{s}$ and
$\alpha$, evenly distributed over $T_{\text{RF}}/T_{2}^{s}\in[6.7,200]$ and
$\alpha\in[0,\pi]$, and numerically optimized $R_{2}^{s,l}$ so that the linear
model (Eq. (23)) closely approximates the $z^{s}$-value at the end of the
pulse as calculated with the generalized Bloch model. This numerical
optimization was performed with the NLsolve.jl package that implements a
trust-region algorithm.23 We approximated
$R_{2}^{s,l}(T_{\text{RF}}/T_{2}^{s},\alpha)$ by cubic B-spline polynomials
with the Interpolations.jl package. Thereafter, we solved Eq. (23) with the
approximated and interpolated $R_{2}^{s,l}(T_{\text{RF}}/T_{2}^{s},\alpha)$ by
taking the matrix exponential, and compared the result during and after RF-
pulses to the solution of the full generalized Bloch-McConnell model (Eqs.
(10)-(13)).

## Results

### Continuous-wave simulations

A generalized-Bloch simulation of $z$-magnetization during a continuous wave
reveals several of the model’s properties. First, after an initial transition,
the magnetization approaches a steady state that is equivalent to the one
predicted by Henkelman’s steady-state model (Eq. (14)) for all three
lineshapes, as visualized in Fig. 2a. This alignment is expected as, for the
case of a steady state, the equivalence of these two models can be shown
mathematically: given a continuous wave ($\omega_{x}=\omega_{1}\sin(\Delta t)$
and $\omega_{y}=\omega_{1}\cos(\Delta t)$, where $\omega_{1}$ is a constant
Rabi frequency and $\Delta$ denotes the off-resonance frequency), and assuming
a steady state ($\partial_{t}z=0$), the generalized Bloch model (Eq. (9))
reduces to

$\begin{split}R_{1}^{s}(1-z_{\text{ss}}^{s})=&\lim_{t\rightarrow+\infty}z_{\text{ss}}^{s}\omega_{1}^{2}\Bigg{(}\cos(\Delta
t)\int_{0}^{t}G(t,\tau)\cos(\Delta\tau)d\tau\\\ &+\sin(\Delta
t)\int_{0}^{t}G(t,\tau)\sin(\Delta\tau)d\tau\Bigg{)}\\\
=&\;\pi\omega_{1}^{2}g(\Delta,T_{2}^{s})z_{\text{ss}}^{s},\end{split}$ (24)

which is equivalent to Henkelman’s steady-state model (Eqs. (14) and (15)).
Here, we used the definition of $G(t,\tau,T_{2}^{s})$ as the Fourier
transformation of $g(\Delta,T_{2}^{s})$.

Figure 2: Simulated dynamics of an isolated semi-solid pool during a
continuous wave. As predicted by theory, we observe good agreement between the
generalized Bloch (Eq. (9)) and Henkelman’s model (Eq. (14)) for the case of a
steady state, as well as between the generalized Bloch model with a Lorentzian
lineshape and the original Bloch model. A comparison to Sled’s theory (Eq.
(21)) reveals deviations between the two models for stronger $\omega_{1}$
fields. The horizontal lines indicate the steady-state predicted by
Henkelman’s theory with respective lineshape. The simulations were performed
with $T_{2}^{s}=10\upmu$s, $R_{1}^{s}=1/$s, and $\Delta/(2\pi)=100$Hz.

Second, Fig. 2 shows good agreement between the generalized Bloch model with a
Lorentzian lineshape and the original Bloch model. This is also anticipated as
the two models are—by virtue of our derivation—mathematically equivalent when
assuming a Lorentzian lineshape. This can be shown by inserting the Green’s
function in Eq. (6) into the generalized Bloch equations (Eqs. (4), (5), and
(9)) and taking the derivative with respect to time. Applying Leibniz’s
integral rule returns the original Bloch equations (Eqs. (1)-(3)).

Third, Fig. 2 compares the generalized Bloch model with Sled’s model for a
Gaussian and a super-Lorentzian lineshape. With an amplitude of
$\omega_{1}/(2\pi)=100$Hz, we observe good agreement between the models (Fig.
2a). For an amplitude of $\omega_{1}/(2\pi)=1000$Hz, however, we find
substantial deviations between the models in the transient phase when assuming
a super-Lorentzian lineshape (Fig. 2b). While Sled’s model (Eq. (21))
describes a decay of the $z^{s}$-magnetization, the generalized Bloch model
(Eq. (9)) captures the angular momentum algebra by modeling the RF-induced
rotation of the magnetization—similar to the Bloch model—which results here in
a negative $z^{s}$-magnetization in the time span $t\gtrapprox 1$ms.

While deviations between Sled’s models and the generalized Bloch model are
present for all lineshapes, they are pronounced the most for the super-
Lorentzian lineshape as the corresponding Green’s function contains more long-
lived components when assuming the same $T_{2}^{s}=10\upmu$s (cf. Fig. 1). A
mathematical comparison of the two models (Eq. (9) vs. Eq. (21)) reveals their
equivalence when assuming $z^{s}(t)=z^{s}(\tau)$ and
$\omega_{1}(\tau)=\omega_{1}(t)$. When this condition is not fulfilled,
however, the subtle difference between Eq. (9) and Eq. (21) can evoke
substantially different spin dynamics, as demonstrated by Fig. 2b.

Figure 3: Simulated $z$-magnetization of an isolated semi-solid spin pool at
the end of an RF-pulse. For a Lorentzian lineshape, we observe good agreement
between the generalized Bloch and the original Bloch model (left column). For
long RF-pulse durations ($T_{\text{RF}}$), we observe good agreement between
all models for all lineshapes. However, for short RF-pulses, only the
generalized Bloch model coincides with the known spin-dynamics, i.e. a
rotation of the magnetization so that
$z^{s}(T_{\text{RF}})/z^{s}(0)=\cos\alpha$, where $\alpha$ denotes the flip
angle. For the super-Lorentzian lineshape, Graham’s single frequency could not
be simulated due to its singularity at $\Delta=0$.

### RF-pulse simulations

The agreement between the generalized Bloch model with a Lorentzian lineshape
and the original Bloch model is confirmed by simulations of the spin dynamics
evoked by RF-pulses (left column of Fig. 3). Fig. 3 further shows good
agreement between all described models for long RF-pulse durations
($T_{\text{RF}}\gg T_{2}^{s}$). However, the figure also demonstrates that all
established models fail to describe the expected rotation induced by short RF-
pulses ($T_{\text{RF}}\ll T_{2}^{s}$). In this limit the pulse approaches a
hard pulse that instantaneously rotates all magnetization and we expect
$z^{s}(T_{\text{RF}})=z^{s}(0)\cos\alpha$. In contrast to established models,
the generalized Bloch model describes such spin dynamics adequately.

The pulse duration at which the generalized Bloch model starts to deviate from
Graham’s and Sled’s models depends on several factors. First, deviations occur
already at much longer $T_{\text{RF}}$ for the super-Lorentzian lineshape
compared to the Lorentzian and Gaussian lineshapes, which can be explained by
the slower signal decay in the corresponding Green’s function (cf. Fig. 1) and
is in line with the observations made for a continuous wave (cf. Fig. 2).
Second, the deviations seem to depend on the flip angle or, more precisely, on
$\omega_{1}$ (cf. the vertical lines in Fig. 3). This analysis indeed reveals
deviations between the models for RF-pulses commonly used in clinical MRI, in
particular for the super-Lorentzian lineshape: e.g., for a $\pi$-pulse with
$T_{\text{RF}}=1$ms, we find a deviation of
$(z_{\text{Sled}}^{s}(T_{\text{RF}})-z_{\text{generalized
Bloch}}^{s}(T_{\text{RF}}))/z^{s}(0)\approx 0.03$. When reducing the pulse
duration to $T_{\text{RF}}=100\upmu$s, which brings the pulse amplitude in the
range of clinical systems’ maximum $B_{1}$-fields, we observe a deviation of
$(z_{\text{Sled}}^{s}(T_{\text{RF}})-z_{\text{generalized
Bloch}}^{s}(T_{\text{RF}}))/z^{s}(0)\approx 0.35$. These simulations were
performed with $T_{2}^{s}=10\upmu$s. The results scale with the value of
$T_{2}^{s}$, i.e. for longer $T_{2}^{s}$ relaxation times, the deviations
become more pronounced.

### Experimental validation

Figure 4: Mono-exponential fits to experimental inversion recovery data of the
MnCl2 and bovine serum albumin (BSA) samples with different inversion pulse
durations $T_{\text{RF}}$ and different inversion times $T_{\text{i}}$. The
data of the MnCl2 sample is well-described by the mono-exponential model
regardless of $T_{\text{RF}}$, as is the data of the BSA sample when inverting
the magnetization with a $22.8\upmu$s pulse. Inverting the magnetization of
the BSA sample with a $912\upmu$s pulse results, however, in a spin-dynamics
that is clearly not mono-exponential.

Figure 5: Inversion recovery NMR experiment of the bovine serum albumin (BSA)
sample with different inversion pulse durations $T_{\text{RF}}$ and different
inversion times $T_{\text{i}}$. This behavior is well-described by the
generalized Bloch model, and better compared to Graham’s spectral model and
Sled’s model. The fitted parameters are shown in Tab. 1.

First, we analyzed the IR curve of each $T_{\text{RF}}$ separately by fitting
a mono-exponential curve (Fig. 4 shows the fits with the shortest and the
longest $T_{\text{RF}}$). The data of the MnCl2 sample is well-described by
the mono-exponential model, regardless of $T_{\text{RF}}$, which is evident by
the small mean relative residual of all fits: $(0.138\pm 0.034)\%$. The
biggest difference between the curves is their inversion efficiency
$-z(0)/z(+\infty)$, which can be explained by $T_{2}^{f}$-relaxation during
the RF-pulse (cf. Fig. 4). When fitting the BSA data with a mono-exponential
model, we do see substantial differences between the individual fits: using an
inversion pulse with $T_{\text{RF}}=22.8\upmu$s, the mono-exponential model
fits the data well (relative residual of 0.13%), which indicates that such
short inversion pulses do indeed invert the magnetization of both spin pools
(free and semi-solid) and we observe an effective $T_{1}$-relaxation of a
fully mixed spin system, i.e. $z^{f}(t)/m_{0}^{f}\approx z^{s}(t)/m_{0}^{s}$
for the entire experiment. In contrast, a pulse with the duration of
$T_{\text{RF}}=912\upmu$s inverts the magnetization of the free pool, but only
saturates the one of the semi-solid pool. As a result, we observe the bi-
exponential decay that is typical for such MT-experiments21 and the mono-
exponential fit is overall poor (relative residual of 2.7%).

Second, we fitted the entire BSA dataset with a global fit using the Bloch-
McConnell model, where the semi-solid pool was described with Graham’s
spectral, Sled’s, or the generalized Bloch model. The fits are depicted in
Fig. 5 and the fitted values are shown in Tab. 1. One can appreciate visually
that the generalized Bloch model fits the measured data much better than
Sled’s and Graham’s models. This is confirmed by the relative residuals: when
fitting the generalized Bloch model, the residual is smaller compared to a
global Bloch fit of the MnCl2 data (0.41% vs. 0.70%). When fitting Graham’s
spectral or Sled’s model to the BSA data, the residual is 4-5 times larger
(1.8% and 2.0%; cf. Tab. 1). Further, the parameters estimated with the
generalized Bloch model (Tab. 1) are in line with literature, while we observe
unrealistically short $T_{2}^{s}$ values for Graham’s and Sled’s model. Yet,
we note that the values are not directly comparable as different studies
measured different samples with different experiments under different
conditions, and analyzed the data with different biophysical models.

Third, we used to the parameters from the generalized Bloch fit to simulate
signals with all three models and we compared the norms of the residuals when
subtracting such simulated signal from the measured one (Fig. 6). This
analysis visualizes that Graham’s and Sled’s models describe the measured
signal particularly poorly for $T_{\text{RF}}\lesssim 300\upmu$s. Yet, the
generalized Bloch model fits the signal better for all inversion pulse
durations measured here, highlighting that even small difference in the models
(e.g., $(z_{\text{Sled}}^{s}(T_{\text{RF}})-z_{\text{generalized
Bloch}}^{s}(T_{\text{RF}}))/z^{s}(0)\approx 0.03$ for a $\pi$-pulse with
$T_{\text{RF}}=1$ms; cf. Fig. 3) can be observed experimentally.

Additional analyses of both samples can be found on
https://jakobasslaender.github.io/MRIgeneralizedBloch.jl/v0.2.11/build_literate/Analyze_NMR_IR_Data.

| Graham’s spectral model | Sled’s model | generalized Bloch model | literature values
---|---|---|---|---
$\frac{||\text{residual}||_{2}}{||\text{signal}||_{2}}(\%)$ | $1.8$ | $2.0$ | $0.41$ |
$m_{0}^{s}$ (%) | $6.99\pm 0.64$ | $7.15\pm 0.38$ | $8.54\pm 0.10$ |  $7.8\pm 0.4$21
$T_{2}^{s}$ ($\upmu$s) | $2.9\pm 5.6$ | $5.0\pm 1.1$ | $12.95\pm 0.61$ |  $10.0\pm 1.0$17
$R_{1}$ (1/s) | $1.179\pm 0.067$ | $1.1782\pm 0.0074$ | $1.1978\pm 0.0015$ |  $1.09\pm 0.03$21
$R_{x}$ (1/s) | $73.9\pm 8.6$ | $60.0\pm 8.1$ | $71.3\pm 2.0$ |  $49.2\pm 5.5$21
Table 1: Parameters of the bovine serum albumin (BSA) sample estimated by
fitting respective model to the data. Selected parameter values reported in
the literature are presented for plausibility. Note, however, that the values
reported here for various models cannot be compared directly to the literature
values, as Ref. 21 used an inversion recovery experiment with
$T_{\text{RF}}=1.5$ms, assuming a Gaussian lineshape and Graham’s single-
frequency model, and Ref. 17 used an off-resonant continuous-wave saturation
experiment. Furthermore, we measured a 15% solution of bovine serum albumin
(BSA), thermally cross-linked, while Ref. 21 chemically cross-linked a 15%
(BSA) solution, and Ref. 17 measured mouse white matter. Lastly, field
strengths and, likely, temperatures differed between the experiments, which
could have influenced $R_{1}$ and $R_{\text{x}}$ further.

Figure 6: Comparison of relative residuals between different models. a:
Residuals of the actual fits, i.e. this analysis uses the biophysical
parameters of respective fit to model the signal. The disadvantage of this
analysis is that residuals at long $T_{\text{RF}}$ are negatively affected by
Graham’s and Sled’s models’ poor description of the signal at short
$T_{\text{RF}}$. b: This problem is overcome by subtracting the measured
signal from signal that is simulated with the estimates from the fit with the
generalized Bloch model. One can observe reduced residuals for Graham’s and
Sled’s models at long $T_{\text{RF}}$. Yet, they are still substantially
larger compared to the ones of the generalized Bloch model throughout, and
most pronounced with an inversion pulse duration of $T_{\text{RF}}\lesssim
300\upmu$s.

### Linear approximation of the generalized Bloch model

The linear approximation achieved approximately a 15,000-fold speed-up. For
example, our current implementation of the full generalized Bloch model took
on average 29ms for a $\pi$-pulse with $T_{\text{RF}}=100\upmu$s on a single
core of a 2015 desktop computer with an Intel Core i5 (I5-6500) processor. In
contrast, the linear approximation, which comprises an evaluation of the
B-spline polynomials and a matrix-exponential, took on average 2.0$\upmu$s. To
give a practical example for the implications: calculating the gradient with
finite differences and using an iterative fitting routine with 30 iterations
for a $256^{3}$ matrix size, would result in a 38 days long computation on 40
cores when using the full generalized Bloch model, while the linear
approximation would take only about 3.8 minutes.

Figure 7: Linearized $T_{2}^{s,l}$ of a super-Lorentzian lineshape. The values
were calculated by solving the generalized Bloch model (Eq. (9)) and comparing
the resulting $z$-magnetization to the linear Bloch-McConnell model (Eq.
(23)).

Key to this speed-up is the pre-computation of the linearized $T_{2}^{s,l}$ as
a function of $T_{2}^{s}$, $T_{\text{RF}}$, and $\alpha$. For a super-
Lorentzian lineshape, the pre-computed values are depicted in Fig. 7. From
this illustration, it becomes apparent that $T_{2}^{s,l}$ depends non-linearly
on all three parameters. Further, we can see that the linearized $T_{2}^{s,l}$
is in the range of $4T_{2}^{s}-8T_{2}^{s}$, where $T_{2}^{s,l}$ is—by
definition of the linear Bloch-McConnell equation—defined by the Lorentzian
lineshape (Eq. (A.25)) and $T_{2}^{s}$ is defined by the super-Lorentzian
lineshape (Eq. (A.27)). This observation is in line with the different decay
characteristics depicted in Fig. 1.

Figure 8: Linear approximation of the generalized Bloch model. a: Plotting the
spin dynamics during a 100$\upmu$s $\pi$-pulse reveals deviations in $z^{s}$
between the generalized Bloch model and its linear approximation. The
linearized relaxation time $T_{2}^{sl}(\alpha,T_{\text{RF}},T_{2}^{s})$ was
chosen such that the approximation matches the full model at the end of the
RF-pulse. The free pool is virtually unaffected by this approximation due to
the comparably slow exchange rate. b: Plotting the spin state at the end of
100$\upmu$s pulses with different flip angles reveals a good agreement between
the generalized Bloch model and its linear approximation.

A comparison of the spin dynamics during an RF-pulse, as modeled by the full
generalized Bloch-McConnell equations (Eqs. (10)-(13)), to the linear
approximation (Eq. (22)) confirms that the free spin-pool is hardly affected
by the approximation, but we do see deviations in the dynamics of the semi-
solid pool, as expected (Fig. 8a). However, by design, the magnetization of
the semi-solid pool simulated with the two models matches at the end of the
RF-pulse. This match also holds true when varying the flip angle, as shown in
Fig. 8b. The normalized root-mean square error of the approximation, measured
at the end of the RF-pulse and averaged over $\alpha\in[0,\pi]$, is $1.5\cdot
10^{-5}$ in the $x^{f}$ component, $6.3\cdot 10^{-6}$ in the $z^{f}$
component, and $2.7\cdot 10^{-4}$ in the $z^{s}$ component.

## Discussion and Conclusions

We have presented a generalization of the Bloch model to non-exponential
decays or, equivalently, to non-Lorentzian lineshapes. We derived and tested
the model in the context of pulsed magnetization transfer1, 2, 10 and
demonstrated that the generalized Bloch model—in contrast to existing pulsed-
MT models10, 11—unifies the original Bloch model and Henkelman’s steady-state
model, while also adequately describing the rotation induced by short RF-
pulses, e.g., the inversion of the magnetization induced by a short
$\pi$-pulse, as observed by Manning et al.12 (cf. Fig. 3).

For small pulse amplitudes ($\omega_{1}\ll 1/T_{2}^{s}$), the generalized
Bloch model shows good agreement with Graham’s and Sled’s models. With an
increasing pulse amplitude, however, we find substantial deviations between
the models (Fig. 3) and NMR experiments confirm that the generalized Bloch
model describes the spin dynamics more accurately (Fig. 5, 6).

Sled’s model (Eq. (21)) and the generalized Bloch model (Eq. (9)) are
equivalent when assuming $z^{s}(t)=z^{s}(\tau)$ and
$\omega_{1}(\tau)=\omega_{1}(t)$. Despite this resemblance, Eqs. (21) and (9)
are mathematically quite different. As Sled’s model has the form
$(\partial_{t}+R_{1})z(t)=-R_{\text{RF}}(t)z(t)+R_{1}$, it describes an
exponential decay with a temporally varying decay rate. The generalized Bloch
model, on the other hand, has the form
$(\partial_{t}+R_{1})z(t)=-\omega_{y}(t)x(t)+R_{1}$, which describes a
rotation in combination with relaxation. This latter feature constitutes the
key property of the generalized Bloch model: explicitly modeling the rotation
of the semi-solid pool’s magnetization ties the magnetization transfer theory
back into the algebra of angular momentum, which forms the basis of the
quantum mechanical description of spin dynamics, as well as the original Bloch
model.

In addition to rotations, a quantum mechanical description entails unitary
evolutions that model the coherence transfer between different spins.
Classical models cannot describe these dynamics to the full extent. Instead,
the generalized Bloch model is based on the assumption that such unitary
transformations average to a net dissipative effect that can be captured by
appropriate Green’s functions. The good agreement of the generalized Bloch
model with the pulsed-MT experiment presented here, but also of Graham’s and
Sled’s models with corresponding experiments that use long pulses with
$\omega_{1}\ll 1/T_{2}^{s}$, gives reason to believe that this assumption is
valid for a heterogeneous mix of large molecules and most clinically feasible
RF-pulses, whose amplitudes commonly vary at a time scale much larger than
$T_{2}^{s}$. A detailed analysis of the generalized Bloch model’s limits of
validity is beyond the scope of this paper, but one would expect that a rapid
series of short pulses with large flip angles could evoke more complex
behavior such as echoes that are not described by the generalized Bloch model.
Applications of the generalized Bloch model beyond magnetization transfer will
also require scrutiny depending on the quantum mechanical processes one aims
to approximate.

Nevertheless, we believe that the generalized Bloch model derived and
validated here is the most complete theory yet formulated for pulsed-MT. Up
until now, the model arguably most commonly used for pulsed-MT is Graham’s
single-frequency approximation. One of the biggest weaknesses of Graham’s
approximation is that the saturation rate is not well-defined for on-resonant
RF-pulses when assuming a super-Lorentzian lineshape (due to its divergence at
the resonant frequency). Gochberg et al. circumvent this problem by
heuristically assuming a Gaussian distribution and relying on the limited
impact of the saturation on their selective inversion recovery MT approach.
21, 24, 25 The bSSFP-based MT approach proposed by Gloor et al. heavily relies
on the saturation of the semi-solid pool. Hence, they do use a super-
Lorentzian lineshape and avoid its singularity by heuristically cutting the
lineshape at 1kHz.26 Like Graham’s spectral model and Sled’s model, the
generalized Bloch model does not face this issue and does not require an
approximation for on-resonant RF-pulses and a super-Lorentzian lineshape. In
addition, the generalized Bloch model provides a more accurate description of
the spin dynamics for short RF-pulses compared to existing models, thus,
extending the well-described experimental design space.

The main disadvantage of the generalized Bloch model is its numerically
challenging structure. With this theory, we foremost aim for a ground truth to
which numerically efficient approximations can be compared. In some cases,
Graham’s and Sled’s models might be sufficiently accurate. For others, we
provided a linear approximation that allows the simulation of the dynamics
during an RF-pulse in about 2$\upmu$s, rendering a broad application of the
generalized Bloch model feasible. The path to establishing a ground truth will
entail further and thorough experimental validation.

##  Appendix A Lineshapes

Throughout this paper we use the following three lineshapes:

$\displaystyle g^{\text{Lorentzian}}(\Delta^{\prime},T_{2}^{s})$
$\displaystyle=\frac{T_{2}^{s}}{\pi}\frac{1}{1+(T_{2}^{s}\Delta^{\prime})^{2}}$
(A.25) $\displaystyle g^{\text{Gaussian}}(\Delta^{\prime},T_{2}^{s})$
$\displaystyle=\frac{T_{2}^{s}}{\sqrt{2\pi}}\exp\left(\frac{-(T_{2}^{s}\Delta^{\prime})^{2}}{2}\right)$
(A.26) $\displaystyle g^{\text{super-L.}}(\Delta^{\prime},T_{2}^{s})$
$\displaystyle=T_{2}^{s}\sqrt{2/\pi}\int_{0}^{1}d\zeta\frac{\exp\left(-2\left(\frac{T_{2}^{s}\Delta^{\prime}}{|3\zeta^{2}-1|}\right)^{2}\right)}{|3\zeta^{2}-1|},$
(A.27)

which are equivalent to the ones reported in Ref. 27.

##  Appendix B Attenuation rates in Graham’s spectral model

Here, we calculate $\langle R_{\text{RF}}\rangle$ for a rectangular RF-pulse
with all three lineshapes defined above. We can define the spectral power
density20 for any RF-pulse $\omega_{1}(t)$ by

$S(\Delta^{\prime},\Delta,\omega_{1}(t))=\frac{1}{\sqrt{2\pi}\;T_{\text{RF}}}\left|\int_{-T_{\text{RF}}/2}^{T_{\text{RF}}/2}\omega_{1}(t)\exp(i(\Delta-\Delta^{\prime})t)dt\right|^{2},$
(A.28)

which describes the spectral distribution of the RF-power. For an on-resonant
rectangular pulse, this Fourier integral returns a squared sinc-function:

$S(\Delta^{\prime},\alpha,T_{\text{RF}})=\frac{\alpha^{2}}{2\pi
T_{\text{RF}}}\text{sinc}^{2}\left(\frac{T_{\text{RF}}\Delta^{\prime}}{2}\right),$
(A.29)

where we assumed an on-resonant RF-pulse ($\Delta=0$) and expressed its
amplitude as a function of its flip angle $\alpha$ and its duration
$T_{\text{RF}}$.

By inserting Eq. (A.29) and the Lorentzian lineshape (Eq. (A.25)) in Eq. (17),
we can analytically solve the integral over $\Delta^{\prime}$ to calculate the
saturation rate of a Lorentzian lineshape for a rectangular RF-pulse:

$R_{\text{RF}}^{\text{Lorentzian}}(\alpha,T_{\text{RF}},T_{2}^{s})=\frac{\alpha^{2}T_{2}^{s}}{T_{\text{RF}}^{2}}\left(\frac{T_{2}^{s}}{T_{\text{RF}}}\exp\left(-\frac{T_{\text{RF}}}{T_{2}^{s}}\right)-\frac{T_{2}^{s}}{T_{\text{RF}}}+1\right).$
(A.30)

For a Gaussian lineshape, this procedure results in

$\begin{split}R_{\text{RF}}^{\text{Gaussian}}&(\alpha,T_{\text{RF}},T_{2}^{s})=\frac{\alpha^{2}T_{2}^{s}}{T_{\text{RF}}^{2}}\\\
&\cdot\Bigg{(}\frac{T_{2}^{s}}{T_{\text{RF}}}\left(\exp\left(-\frac{T_{\text{RF}}^{2}}{2{T_{2}^{s}}^{2}}\right)-1\right)+\sqrt{\frac{\pi}{2}}\text{erf}\left(\frac{T_{\text{RF}}}{\sqrt{2}T_{2}^{s}}\right)\Bigg{)}.\end{split}$
(A.31)

Here, $\text{erf}(x)$ denotes the error function. For a super-Lorentzian
lineshape, Fubini’s theorem tells us that we can switch the order of the two
integrals. Solving the integral over $\Delta^{\prime}$ results in

$R_{\text{RF}}^{\text{Super-L.}}(\alpha,T_{\text{RF}},T_{2}^{s})=\frac{\alpha^{2}T_{2}^{s}}{T_{\text{RF}}^{2}}f(T_{\text{RF}}/T_{2}^{s})$
(A.32)

with

$\begin{split}f(\tau)&=\int_{0}^{1}\frac{\sqrt{2\pi}}{|3\zeta^{2}-1|}\text{erf}\left(\frac{\tau}{\sqrt{8}}|3\zeta^{2}-1|\right)\\\
&+\frac{4}{\tau|3\zeta^{2}-1|^{2}}\left(\exp\left(-\frac{\tau^{2}|3\zeta^{2}-1|^{2}}{8}\right)-1\right)d\zeta.\end{split}$
(A.33)

For all three lineshapes, the saturation rate has the form
$\alpha^{2}T_{2}^{s}/T_{\text{RF}}^{2}\cdot f(T_{\text{RF}}/T_{2}^{s})$, where
the difference between the lineshapes lies in different functions $f(\tau)$.

##  Appendix C Spectral analysis of Graham’s single frequency approximation

In order to analyze the limitations of Graham’s single frequency
approximation, we plotted in Fig. A.9 a Lorentzian, a Gaussian, and a super-
Lorentzian lineshape along with the spectral power density functions of
rectangular RF-pulses. For a Lorentzian or Gaussian lineshape in combination
with an on-resonant, rectangular RF-pulse with $T_{\text{RF}}=500\upmu$s,
Graham’s approximations appear justified: neither of those two lineshapes
exhibit substantial changes throughout the spectral profile of the RF-pulse.
We can, hence, approximate the lineshape by a constant value. Reducing
$T_{\text{RF}}$ to $100\upmu$s pushes the limits of such an approximation: the
lineshapes show some variations over the larger spectral profile and we can
expect some errors when using this approximation. For a super-Lorentzian
lineshape, this approximation cannot be made at all: due to the singularity on
resonance ($\Delta^{\prime}=0$), the lineshape varies throughout the power
spectral density, regardless of the pulse shape and duration.

Figure A.9: Spectral analysis of the three lineshapes discussed in this paper
and RF-pulses with different pulse durations $T_{\text{RF}}$. The Lorentzian
and Gaussian lineshapes $g(\Delta^{\prime},T_{2}^{s})$ are constant throughout
the spectral power density $S(\Delta^{\prime},\alpha,T_{\text{RF}})$ of an on-
resonant rectangular pulse with $T_{\text{RF}}=500\upmu$s and, to a lesser
degree, with $T_{\text{RF}}=100\upmu$s. This assumption is not valid for a
super-Lorentzian lineshape due to the divergence of the lineshape on-
resonance. All lineshapes were simulated with $T_{2}^{s}=10\upmu$s and we set
$\alpha=\pi$ for both RF-pulses.

##  Appendix D Acknowledgements

We would like to thank Rafael P. Bruschweiler for fruitful discussions.

##  Appendix E Data availability statement

In order to ensure reproducibility, we provide the following resources:

  * •

Highly efficient implementations of the generalized Bloch model and its linear
approximation can be found on
https://github.com/JakobAsslaender/MRIgeneralizedBloch.jl. They are written in
the open source language Julia and we registered the package
"MRIgeneralizedBloch.jl," which facilitates the installation via Julia’s
package manager.

  * •

Scripts to reproduce all simulations, data analyses, and figures can be found
in the documentation of the package:
https://JakobAsslaender.github.io/MRIgeneralizedBloch.jl/v0.2.11. They render
the code in HTML format with interactive figures and link to Jupyter notebooks
that can be launched in binder, allowing for the replication and modification
of all simulations in a browser without any local installations.

  * •

All NMR data are stored in a separate github repository
https://github.com/JakobAsslaender/MRIgeneralizedBloch_NMRData and are loaded
and analyzed by above mentioned scripts.

All figures and tables in this paper were created with the version 0.2.11
(commit hash b185969) of the MRIgeneralizedBloch.jl package and match the
version’s documentation. The commit hash of the data used for the preparation
of this paper is 77950c4.

## References

  * 1 Wolff SD, Balaban RS. Magnetization transfer contrast (MTC) and tissue water proton relaxation in vivo. Magn. Reson. Med. 1989; 10(1): 135–144.  doi: 10.1002/mrm.1910100113
  * 2 Henkelman RM, Huang X, Xiang QSS, Stanisz GJ, Swanson SD, Bronskill MJ. Quantitative Interpretation of Magnetization Transfer. Magn. Reson. Med. 1993; 29(6): 759–766.
  * 3 Melki PS, Mulkern RV. Magnetization transfer effects in multislice RARE sequences. Magn. Reson. Med. 1992; 24(1): 189–195.  doi: 10.1002/mrm.1910240122
  * 4 Bieri O, Scheffler K. On the Origin of Apparent Low Tissue Signals in Balanced SSFP.. Magn. Reson. Med. 2006; 56(5): 1067–74.  doi: 10.1002/mrm.21056
  * 5 McConnell HM. Reaction Rates by Nuclear Magnetic Resonance. J. Chem. Phys. 1958; 28(3): 430–431.  doi: 10.1063/1.1744152
  * 6 Morrison C, Mark Henkelman R. A Model for Magnetization Transfer in Tissues. Magn. Reson. Med. 1995; 33(4): 475–482.  doi: 10.1002/mrm.1910330404
  * 7 Bloembergen N, Purcell EM, Pound RV. Relaxation Effects in Nuclear Magnetic Resonance Absorption. Phys. Rev. 1948; 73(7): 679–712.  doi: 10.1103/PhysRev.73.679
  * 8 Bloch F. Nuclear induction. Phys. Rev. 1946; 70(7-8): 460–474.  doi: 10.1103/PhysRev.70.460
  * 9 Hahn EL. Spin Echoes. Phys. Rev. 1950; 80(4): 580–594.  doi: 10.1103/PhysRev.80.580
  * 10 Graham SJ, Henkelman RM. Understanding pulsed magnetization transfer. J. Magn. Reson. Imaging 1997; 7(5): 903–912.  doi: 10.1002/jmri.1880070520
  * 11 Sled JG, Pike GB. Quantitative Interpretation of Magnetization Transfer in Spoiled Gradient Echo MRI Sequences. J. Magn. Reson. 2000; 145(1): 24–36.  doi: 10.1006/jmre.2000.2059
  * 12 Manning AP, MacKay AL, Michal CA. Understanding aqueous and non-aqueous proton T1 relaxation in brain. J. Magn. Reson. 2021; 323: 106909\.  doi: 10.1016/j.jmr.2020.106909
  * 13 Rabi I, Zacharias J, Millman S, Kusch P. A New Method of Measuring Nuclear magnetic Moment. Phys. Rev. 1938; 53(February): 318.  doi: 10.1103/PhysRev.53.318
  * 14 Ernst RR, Anderson WA. Application of Fourier Transform Spectroscopy to Magnetic Resonance. Rev. Sci. Instrum. 1966; 37(1): 93–102.  doi: 10.1063/1.1719961
  * 15 Ernst R, Aue W, Bartholdi E, Höhener A, Schäublin S. Equivalence of Fourier Spectroscopty and Slow Passage in Nuclear Magnetic Resonance. In: No. iii. Elsevier. 1974 (pp. 47–60)
  * 16 Wennerström H. Proton nuclear magnetic resonance lineshapes in lamellar liquid crystals. Chem. Phys. Lett. 1973; 18(1): 41–44.  doi: 10.1016/0009-2614(73)80333-1
  * 17 Stanisz GJ, Odrobina EE, Pun J, et al. T1, T2 relaxation and magnetization transfer in tissue at 3T.. Magn Reson Med 2005; 54(3): 507–512.  doi: 10.1002/mrm.20605
  * 18 Rackauckas C, Nie Q. DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia. J. Open Res. Softw. 2017; 5.  doi: 10.5334/jors.151
  * 19 Assländer J, Novikov DS, Lattanzi R, Sodickson DK, Cloos MA. Hybrid-state free precession in nuclear magnetic resonance. Nat. Commun. Phys. 2019; 2(1): 73.  doi: 10.1038/s42005-019-0174-0
  * 20 Stremler FG. Introduction to communication systems.  Addison Wesley Publishing Company .  1982\.
  * 21 Gochberg DF, Gore JC. Quantitative imaging of magnetization transfer using an inversion recovery sequence. Magn. Reson. Med. 2003; 49(3): 501–505.  doi: 10.1002/mrm.10386
  * 22 Wimperis S. Broadband and narrowband composite excitation sequences. J. Magn. Reson. 1990; 86(1): 46–59.  doi: 10.1016/0022-2364(90)90210-Z
  * 23 Nocedal J, Wright S. Numerical optimization.  Springer Science & Business Media .  2006\.
  * 24 Gochberg DF, Gore JC. Quantitative magnetization transfer imaging via selective inversion recovery with short repetition times. Magn. Reson. Med. 2007; 57(2): 437–441.  doi: 10.1002/mrm.21143
  * 25 Dortch RD, Li K, Gochberg DF, et al. Quantitative magnetization transfer imaging in human brain at 3 T via selective inversion recovery. Magn. Reson. Med. 2011; 66(5): 1346–1352.  doi: 10.1002/mrm.22928
  * 26 Gloor M, Scheffler K, Bieri O. Quantitative magnetization transfer imaging using balanced SSFP.. Magn. Reson. Med. 2008; 60(3): 691–700.  doi: 10.1002/mrm.21705
  * 27 Morrison C, Stanisz G, Henkelman R. Modeling Magnetization Transfer for Biological-like Systems Using a Semi-solid Pool with a Super-Lorentzian Lineshape and Dipolar Reservoir. J. Magn. Reson. Ser. B 1995; 108(2): 103–113.  doi: 10.1006/jmrb.1995.1111


# Diophantine $D(n)$-quadruples in $\mathbb{Z}[\sqrt{4k+2}]$

Kalyan Chakraborty, Shubham Gupta and Azizul Hoque  K. Chakraborty @Kerala
School of Mathematics, Kozhikode-673571, Kerala, India<EMAIL_ADDRESS>S. Gupta @Harish-Chandra Research Institute, A CI of Homi Bhabha National
Institute, Chhatnag Road, Jhunsi, Prayagraj - 211019<EMAIL_ADDRESS>A. Hoque @Department of Mathematics, Rangapara College, Rangapara,
Sonitpur-784505, Assam, India<EMAIL_ADDRESS>
###### Abstract.

Let $\mathbb{Z}[\sqrt{d}]$ be a quadratic ring of integers. For a given
$n\in\mathbb{Z}[\sqrt{d}]$, a set of $m$ non-zero distinct elements in
$\mathbb{Z}[\sqrt{d}]$ is called a Diophantine $D(n)$-$m$-tuple in
$\mathbb{Z}[\sqrt{d}]$ if product of any two of them plus $n$ is a square in
$\mathbb{Z}[\sqrt{d}]$. Assume that $d\equiv 2\pmod{4}$ such that the
equations $x^{2}-dy^{2}=-1$ and $x^{2}-dy^{2}=6$ are solvable in integers. In
this paper, we prove that there exist infinitely many Diophantine
$D(n)$-quadruples in $\mathbb{Z}[\sqrt{d}]$ for $n=4m+4k\sqrt{d}$ with
$m,k\in\mathbb{Z}$ satisfying $(m,k)\not\equiv(5,3),(11,3)\pmod{(12,6)}$.
Also, we prove a similar result for $n=(4m+2)+4k\sqrt{d}$ when
$(m,k)\not\equiv(9,3),(0,0)\pmod{(12,6)}$. At the end, we provide some
examples supporting the existence of Diophantine quadruples in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$ for the above exceptions.

###### Key words and phrases:

Diophantine quadruples; Pellian equations; Quadratic fields

###### 2020 Mathematics Subject Classification:

11D09; 11R11

##  1\. Introduction

A set $\\{a_{1},a_{2},\cdots,a_{m}\\}$ of $m$ distinct positive integers is
called a Diophantine $m$-tuple with the property $D(n)$ (or simply
$D(n)$-$m$-tuple) for a given non-zero integer $n$, if $a_{i}a_{j}+n$ is a
perfect square for all $1\leq i<j\leq m$. For $n=1$, such a $m$-tuple is
called Diophantine $m$-tuple instead of Diophantine $m$-tuple with the
property $D(1)$. The problem of constructing such tuples was first studied by
Diophantus of Alexandria, who found a quadruple of rational numbers, viz.
$\\{1/16,33/16,17/4,105/16\\}$ with the property ‘the product of any two of
them increased by unity is the square of a rational number’. However, it was
Fermat who first found a Diophantine quadruple, viz. $\\{1,3,8,120\\}$. Later,
Baker and Davenport [3] proved that Fermat’s quadruple can not be extended to
Diophantine quintuple. Dujella [10] proved in 2004 that there is no
Diophantine sextuple and that there are at most finitely many Diophantine
quintuples. Recently, He, Togbé and Ziegler [15] proved that there is no
Diophantine quintuple, and thus they solved a long-standing open problem. On
the other hand, Bonciocat, Cipu and Mignotte [16] proved the conjecture of
Dujella [8], which states that there is no Diophantine $D(-1)$-quadruple. It
is also known due to Bliznac Trebješanin and Filipin [4] that there is no
Diophantine $D(4)$-quintuple. A brief survey on this topic can be found in
[12]. We also refer [5, 7, 11, 13] to the reader for more information about
Diophantine $D(n)$-$m$-tuples.

Let us look at the above concept in a general set up. Let $\mathcal{R}$ be a
commutative ring with unity. For a given $n\in\mathcal{R}$, a set
$\\{a_{1},a_{2},\cdots,a_{m}\\}\subset\mathcal{R}\setminus\\{0\\}$ is called a
Diophantine $m$-tuple in $\mathcal{R}$ with the property $D(n)$ (or simply
simply $D(n)$-$m$-tuple in $\mathcal{R}$) if $a_{i}a_{j}+n$ is a perfect
square in $\mathcal{R}$ for all $1\leq i<j\leq m$. For the ring of integers
$\mathcal{O}_{K}$ of any imaginary quadratic fields, Adžaga [1] proved that
there is no Diophantine $m$-tuples with the property $D(1)$ in
$\mathcal{O}_{K}$ when $m\geq 42$, while Gupta [14] proved that there is no
such tuple with the property $D(-1)$ when $m\geq 37$. On the other hand, it is
more interesting to see that Diophantine $D(n)$-quadruples are related to the
representations of $n$ by the binary quadratic form $x^{2}-y^{2}$. In
particular, Dujella [8] proved that a Diophantine quadruple in $\mathbb{Z}$
with the property $D(n)$ exists if and only if $n$ can be written as the
difference of two squares, up to finitely many exceptions. Later, Dujella [9]
proved the above fact in the case of Gaussian integers. Further, the above
property also holds for the ring of integers of $\mathbb{Q}(\sqrt{d})$ for
certain $d\in\mathbb{Z}$ (see, [2, 17, 18, 19, 22, 24]). These results
motivate Franušić and Jadrijević to post the following conjecture.

######  Conjecture 1.1 ([21, Conjecture 1]).

Let $\mathcal{R}$ be a commutative ring with unity $1$ and
$n\in\mathcal{R}\setminus\\{0\\}$. Then a quadruple with the property $D(n)$
exists if and only if $n$ can be written as a difference of two squares in
$\mathcal{R}$, up to finitely many exceptions of $n$.

This conjecture was verified for rings of integers of certain number fields
(cf. [2, 17, 18, 19, 20, 21, 22, 23, 24]). Very recently, the authors gave
some counter-examples to this conjecture in [6]. Precisely, they showed that
Conjecture 1.1 does not hold for $\mathbb{Z}[\sqrt{10}]$.

We now fix some notations to proceed further.

  * •

$(a,b)=a+b\sqrt{d}$,

  * •

Let $\alpha=(a,b)$. We define norm of this element $\alpha$ by
$Nm(\alpha):=(a,b)(a,-b),$

  * •

$(x,y)\equiv(a,b)\text{~{}(mod~{}}(c,d))$ means that $x\equiv
a\pmod{c}\text{~{}~{}~{}and~{}~{}~{}}y\equiv b\pmod{d}.$

Assume that $d\equiv 2\pmod{4}$ is a square-free integer and
$m,k\in\mathbb{Z}$. We define the sets $\mathcal{S}$ and $\mathcal{T}$ as
follows:

$\displaystyle\mathcal{S}:=$
$\displaystyle\\{(4m,4k+1),(4m,4k+2),(4m,4k+3),(4m+1,4k+1),(4m+1,4k+3),(4m+2,$
$\displaystyle 4k+1),(4m+2,4k+3),(4m+3,4k+1),(4m+3,4k+3)\\},$
$\displaystyle\mathcal{T}:=$
$\displaystyle\\{(4m,4k),(4m+1,4k),(4m+1,4k+2),(4m+2,4k),(4m+2,4k+2),(4m+3,4k),$
$\displaystyle(4m+3,4k+2)\\}.$

It is easy to check that if $n\in\mathbb{Z}[\sqrt{d}]$ then
$n\in\mathcal{S}\cup\mathcal{T}$. In [17], Franušić proved the following
theorem.

######  Theorem A.

Let $d\equiv 2\pmod{4}$ be a square-free integer. Then for any
$n\in\mathcal{S}$, there does not exist any Diophantine $D(n)$-quadruple in
$\mathbb{Z}[\sqrt{d}]$.

Thus, it would be interesting to ask ‘whether there exists any Diophantine
quadruple in $\mathbb{Z}[\sqrt{d}]$ for $n\in\mathcal{T}$’. In [6], the
authors proved that there are infinitely many Diophantine $D(n)$-quadruple in
$\mathbb{Z}[\sqrt{d}]$ for $n\in\mathcal{T}\setminus\\{(4m,4k),(4m+2,4k)\\}$
with $m,k\in\mathbb{Z}$. In this paper, we consider the above problem for the
remaining two values of $n$. In other words, we prove there are infinitely
many Diophantine $D(n)$-quadruple in $\mathbb{Z}[\sqrt{d}]$ for
$n\in\\{(4m,4k),(4m+2,4k)\\}$ with $m,k\in\mathbb{Z}$. Here $d\equiv
2\pmod{4}$ is a square-free integer with the property that the equations

$x^{2}-dy^{2}=-1$ (1.1)

and

$x^{2}-dy^{2}=6$ (1.2)

are solvable in integers. More precisely, we prove:

######  Theorem 1.1.

Assume that $d\equiv 2\pmod{4}$ is a square-free positive integer and the
equations (1.1) and (1.2) are solvable in integers. Let $n=(4m,4k)$ with
$m,k\in\mathbb{Z}$ such that $(m,k)\not\equiv(5,3),(11,3)\pmod{(12,6)}$. Then
there exist infinitely many Diophantine $D(n)$-quadruple in
$\mathbb{Z}[\sqrt{d}]$.

######  Theorem 1.2.

Let $d$ be as in Theorem 1.1. Then for $n=(4m+2,4k)$ with $m,k\in\mathbb{Z}$,
there exist infinitely many Diophantine $D(n)$-quadruple in
$\mathbb{Z}[\sqrt{d}]$, except for $(m,k)\in\\{(12r+9,6s+3),(12r,6s)\\}$ with
$r,s\in\mathbb{Z}$.

##  2\. Preliminaries

We begin with the following lemma which follows from the definition of
Diophantine $D(n)$-quadruples in $\mathbb{Z}[\sqrt{d}]$.

######  Lemma 2.1.

For a square-free $d\in\mathbb{Z}$, let $\\{a_{1},a_{2},a_{3},a_{4}\\}$ be a
quadruple with the property $D(n)$. Then for any non-zero
$w\in\mathbb{Z}[\sqrt{d}]$, $\\{wa_{1},wa_{2},wa_{3},wa_{4}\\}$ is a quadruple
in $\mathbb{Z}[\sqrt{d}]$ with the property $D(w^{2}n)$.

The next lemma tells helps to find the conditions under which the set
$\\{a,b,a+b+2r,a+4b+4r\\}$ forms a Diophantine $D(n)$-quadruple in
$\mathbb{Z}[\sqrt{d}]$ for any $n\in\mathbb{Z}[\sqrt{d}]$.

######  Lemma 2.2.

[6, Lemma 2.5] Let $n\in\mathbb{Z}[\sqrt{d}]$. A set
$\\{a,b,a+b+2r,a+4b+4r\\}$ having all non-zero and distinct elements form a
Diophantine 4-tuple in $\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$ if
$ab+n=r^{2}$ and $3n=\alpha_{1}\alpha_{2}$, where $\alpha_{1}=a+2r+\alpha$ and
$\alpha_{2}=a+2r-\alpha$ for some $a,b,r,\alpha\in\mathbb{Z}[\sqrt{d}]$.

We also need the next two lemmas in order to apply Lemma 2.2 in the proofs of
Theorems 1.1 and 1.2. Lemma 2.3 is useful for the factorization of $3n$ in
$\mathbb{Z}[\sqrt{d}]$, while Lemma 2.4 is useful for verifying the non-zero
and distinct elements.

######  Lemma 2.3 ([6, Lemma 3.1]).

Let $d\equiv 2\pmod{4}$ is a square-free integer such that the equations (1.1)
and (1.2) are solvable in integers. Then for the ring of integers
$\mathbb{Z}[\sqrt{d}]$, the following statements hold:

  * (i)

All the elements of norm $1$ have the form $(6a_{1}\pm 1,6b_{1})$ and there
are infinitely many such elements;

  * (ii)

All the elements of norm $-1$ have the form $(6a\pm 3,6b\pm 1)$ and there are
infinitely many such elements;

  * (iii)

$d\equiv 10\pmod{48}$;

  * (iv)

All the elements of norm $6$ have the form $(12M\pm 4,6N\pm 1)$ and there are
infinitely many such elements;

  * (v)

All the elements of norm $-6$ have the form $(12M\pm 2,6N\pm 1)$ and there are
infinitely many such elements; where $a,b,a_{1},b_{1},M,N\in\mathbb{Z}$.

######  Lemma 2.4 ([6, Lemma 2.4]).

Assume that
$a_{1},a_{2},b_{1},b_{2},c_{1},c_{2},d_{1},d_{2},e_{1}\in\mathbb{Z}$ with
$a_{1}a_{2}b_{1}\neq 0$. Then the following system of simultaneous equations

$\displaystyle\begin{cases}a_{1}x^{2}+b_{1}y^{2}+c_{1}x+d_{1}y+e_{1}=0,\\\
a_{2}xy+b_{2}x+c_{2}y+d_{2}=0\end{cases}$ (2.1)

has only finitely solutions in integers.

##  3\. Proof of Theorem 1.1

We begin this proof with $n=(4m,4k)$. Here, $3n=3(4m,4k)=6(2m,2k)$ and we
choose $\alpha_{1}=6$ and $\alpha_{2}=(2m,2k)$. Then utilizing Lemma 2.2, we
get

$a+2r=(m+3,k).$ (3.1)

Case I: $m,k$ are odd.
Let us consider $a=(6a_{1}+1,6b_{1})$ with $Nm(a)=1$, by Lemma 2.3 $(i)$.
There are infinitely many choices for such values of $a$ by Lemma 2.3 $(i)$.
Thus (3.1) implies that

$2r=(m+2-6a_{1},k-6b_{1}).$

In this case, both $m$ and $k$ are even, then $r\in\mathbb{Z}[\sqrt{d}]$, and
hence

$r=(m/2+1-3a_{1},k/2-3b_{1}).$

Using these $a$ and $r$, we get $b$ from the equation $ab+n=r^{2}$, which is

$b=((m/2+1-3a_{1})^{2}+d(k/2-3b_{1})^{2}-4m,2(m/2+1-3a_{1})(k/2-3b_{1})-4k)(6a_{1}+1,-6b_{1}).$

Since both $m,k$ are even, so the above choices of $a,b$ and $r$ provide
infinitely many quadruples $\\{a,b,a+b+2r,a+4b+4r\\}$ in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(4m,4k)$, by using Lemma 2.2. The
condition of non-zero distinct elements can be checked using Lemma 2.4.
Analogusly, for other cases of $m$ and $k$, Diophantine 4-tuple in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$ can be obtained.
Case II: $m$ is odd and $k$ is even.
In this case, we choose $a=2(6a_{1}+1,6b_{1})$ with $Nm((6a_{1}+1,6b_{1}))=1$
and then from (3.1) we get

$2r=(m+1-12a_{1},k-12b_{1}).$

Since $m$ is odd and $k$ is even, so $m=2m_{1}+1$ and $k=2k_{1}$, for some
$m_{1},k_{1}\in\mathbb{Z}$. Thus, we have

$r=(m_{1}+1-6a_{1},k_{1}-6b_{1}).$

We utilize $ab+n=r^{2}$ to get

$b=\dfrac{1}{2}\left((m_{1}+1-6a_{1})^{2}-d(k_{1}-6b_{1})^{2}-4m,2(m_{1}+1-6a_{1})(k_{1}-6b_{1})-4k\right)(6a_{1}+1,-6b_{1}).$

For $b\in\mathbb{Z}[\sqrt{d}]$, $m_{1}$ must be chosen as odd, and thus these
choices of $a,b$ and $r$ provide infinitely many quadruples
$\\{a,b,a+b+2r,a+4b+4r\\}$ in $\mathbb{Z}[\sqrt{d}]$ with the property
$D(4m,4k)$, by using Lemma 2.2.

Now for even $m_{1}$, we choose $a=4(6a_{1}+1,6b_{1})$ with
$Nm((6a_{1}+1,6b_{1}))=1$, then (3.1) gives

$r=(m_{1}-12a_{1},k_{1}-12b_{1}).$

As before by $ab+n=r^{2}$, we get

$b=\dfrac{1}{4}\left((m_{1}-12a_{1})^{2}+d(k_{1}-12b_{1})^{2}-4m,2(m_{1}-12a_{1})(k_{1}-12b_{1})-4k\right)(6a_{1}+1,-6b_{1}).$

For $b\in\mathbb{Z}[\sqrt{d}]$, then $k_{1}$ must also be even. One can check
that except $a+4b+4r=0$ and $a+2r=0$, the condition of non-zero and distinct
elements of the set which is given in Lemma 2.2 can be verified by Lemma 2.4.
Let $a+4b+4r=0$, then we must have $n=(1,0)$ or $(36,0)$. Also if $a+2r=0$,
then $n=-12$. However, we have infinitely many quadruples with the property
$D(n)$ for $n=(1,0),(36,0).$ Therefore using these values of $a,b$ and $r$, we
can generate infinitely many quadruples with the property $D(n)$ by Lemma 2.2,
except for $n=-12$.

Let $n=-12.$ We factorize $3n$ in the following way

$3n=3\times(-12)=-36=-18\times 2.$

Now we take $\alpha_{1}=-18$ and $\alpha_{2}=2$. From Lemma 2.2,

$a+2r=-8.$

We choose $a=4(6a_{1}+1,6b_{1})$ with $Nm((6a_{1}+1,6b_{1}))=1$, and then we
get

$r=(-6-12a_{1},-12b_{1}).$

Using $ab+n=r^{2}$,

$b=\dfrac{((-6-12a_{1},-12b_{1})^{2}+12)(6a_{1}+1,-6b_{1})}{4}.$

After simplifying the above equation, we get

$b=((3+6a_{1},6b_{1})^{2}+3)(6a_{1}+1,-6b_{1}).$

Since, due to existence of infintely many $a$’s (by Lemma 2.3 $(i)$), there
exist infinitely many quadruples with the property $D(-12)$, by using Lemma
2.2. Here the condition of non-zero and distinct elements can be checked by
Lemma 2.4.

We now consider another possibility of $m_{1}$ and $k_{1}$ which are even and
odd, respectively then $n$ will turn out to the form
$(16m+4,16k+8)=2^{2}(4m+1,4k+2)$.

As we have a quadruple with the property $D(4m+1,4k+2)$ (see [[6] Theorem
1.1]), so that by Lemma 2.1 one gets a quadruple with the property
$D(16m+4,16k+8)$. This concludes that there are infinitely many quadruples in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(4m,4k),$ for $m$ is odd and $k$ is
even.
Case III: $m$ is even and $k$ is odd.
For this case, we choose:

$a=(6a_{1}+3,6b_{1}+1)$

with $Nm(a)=-1$, by Lemma 2.3 $(ii)$. Then using (3.1) we get

$r=(m/2-3a_{1},(k-1)/2-3b_{1}).$

As before utilizing $ab+n=r^{2}$, we get

$b=((m/2-3a_{1})^{2}+d((k-1)/2-3b_{1})^{2}-4m,2(m/2-3a_{1})((k-1)/2-3b_{1})-4k)(-6a_{1}-3,6b_{1}+1).$

Case IV: Both $m$ and $k$ are odd.
Since $m$ and $k$ are odd, so $n$ is of the form $(8m+4,8k+4)$. Then

$3n=6(4m+2,4k+2).$

Here we choose $\alpha_{1}=6$ and $\alpha_{2}=(4m+2,4k+2)$, which imply that

$a+2r=(2m+4,2k+1).$ (3.2)

Applying Lemma 2.3 $(iv)$, we set $a=(12M+4,6N+1)$ with $Nm(a)=6$ and thus
(3.2) implies that

$r=(m-6M,k-3N).$

Employing $ab+n=r^{2}$ and $d\equiv 10\pmod{48}$ (by Lemma 2.3 $(iii)$), we
get

$b=\frac{1}{6}\left((m-6M)^{2}+d(k-3N)^{2}-8m-4,2(m-6M)(k-3N)-8k-4\right)(12M+4,-6N-1).$

For existence of $b$ in $\mathbb{Z}[\sqrt{d}]$, we must have

$(m,k)\equiv(0,1),(2,2),(0,0),(4,2),(4,1)\pmod{(6,3)}.$

Again, we choose $a=(12M+4,6N-1)$ such that $Nm(a)=6$. Then using (3.2), we
get

$r=(m-6M,k-3N+1).$

As before utilizing $ab+n=r^{2}$, one gets

$\displaystyle b=$
$\displaystyle\dfrac{((m-6M)^{2}+d(k-3N+1)^{2}-8m-4,2(m-6M)(k-3N+1)-8k-4)}{6}$
$\displaystyle\times(12M+4,-6N+1).$

We require $b$ in $\mathbb{Z}[\sqrt{d}]$, so that

$(m,k)\equiv(2,0),(4,0),(0,2)\pmod{(6,3)}.$

Analogously we choose $a=(12M+2,6N+1)$ with $Nm(a)=-6$, by Lemma 2.3 $(v)$. In
this case, (3.2) implies

$r=(m+1-6M,k-3N),$

and we use $ab+n=r^{2}$ to get,

$\displaystyle b=$
$\displaystyle\dfrac{((m+1-6M)^{2}+d(6N+1)^{2}-8m-4,2(m+1-6M)(6N+1)-8k-4)}{-6}$
$\displaystyle\times(12M+2,-6N-1).$

As before, $b$ to be in $\mathbb{Z}[\sqrt{d}]$, one gets

$(m,k)\equiv(1,0),(3,2),(1,1),(5,2)\pmod{(6,3)}.$

We now choose $a=(12M+2,6N-1)$ satisfying $Nm(a)=-6$, by Lemma 2.3 $(v)$ and
then applying (3.2), we get

$r=(m-6M+1,k-3N+1),$

and thus $ab+n=r^{2}$ gives

$\displaystyle b=$
$\displaystyle\dfrac{((m-6M+1)^{2}+d(k-3N+1)^{2}-8m-4,2(m-6M+1)(k-3N+1)-8k-4)}{-6}$
$\displaystyle\times(12M+2,-6N+1).$

For $b\in\mathbb{Z}[\sqrt{d}]$, one must choose

$(m,k)\equiv(3,0),(1,2)\pmod{(6,3)}.$

Again choose $a=(12M-2,6N+1)$ with $Nm(a)=-6$, then (3.2) implies

$r=(m+3-6M,k-3N)$

and

$\displaystyle b=$
$\displaystyle\dfrac{((m+3-6M)^{2}+d(k-3N)^{2}-8m-4,2(m+3-6M)(k-3N)-8k-4)}{-6}$
$\displaystyle\times(12M-2,-6N-1).$

Now for $b\in\mathbb{Z}[\sqrt{d}]$, we get

$(m,k)\equiv(3,1)\pmod{(6,3)}.$

Let $a=(12M-2,6N-1)$ satisfying $Nm(a)=-6$. Then (3.2) gives us

$r=(m+3-6M,k-3N+1)$

and hence

$\displaystyle b=$
$\displaystyle\dfrac{((m+3-6M)^{2}+d(k-3N+1)^{2}-8m-4,2(m+3-6M)(k-3N+1)-8k-4)}{-6}$
$\displaystyle\times(12M-2,-6N+1).$

Now if

$(m,k)\equiv(5,0)\pmod{(6,3)},$

then $b\in\mathbb{Z}[\sqrt{d}]$.

##  4\. Proof of Theorem 1.2

The proof is very similar to proof of the previous theorem. We divide the
proof into some parts depending on the values of $n=(4m+2,4k)$.
Case I:  $2\mid m$ and $2\nmid k$.

We factorize $3n$ as follows:

$3n=6(2m+1,2k)=(12M+4,-6N-1)(12M+4,6N+1)(2m+1,2k).$

Assume that

$\begin{cases}\alpha_{1}=(12M+4,-6N-1),\\\
\alpha_{2}=(24Mm+12M+8m+4+d(12Nk+2k),24Mk+8k+12Nm+2m+6N+1).\end{cases}$

Then by Lemma 2.2, one gets

$a+2r=(12Mm+12M+2+4m+2+d(6Nk+k),12Mk+4k+6Nm+m).$ (4.1)

We first take $a=2(6a_{1}+1,6b_{1})$ such that $Nm((6a_{1}+1,6b_{1}))=1$, by
Lemma 2.3 $(i)$ then

$r=(6Mm+6M+2m+1+(d/2)(6Nk+k)-6a_{1},6Mk+2k+3Nm+(m/2)-6b_{1}).$

We now use $ab+n=r^{2}$ to get

$\displaystyle b$
$\displaystyle=\dfrac{1}{2}\Big{\\{}((6Mm+6M+2m+1+(d/2)(6Nk+k)-6a_{1})^{2}+d(6Mk+2k+3Nm+(m/2)$
$\displaystyle-6b_{1})^{2}-4m-2,2(6Mm+6M+2m+1+(d/2)(6Nk+k)-6a_{1})(6Mk+2k+3Nm$
$\displaystyle+(m/2)-6b_{1})-4k)\times(6a_{1}+1,-6b_{1})\Big{\\}}.$

Since $2\mid m$ and $2\nmid k$, so that $r,b\in\mathbb{Z}[\sqrt{d}]$. Thus
Lemma 2.2 ensures that these values of $a,b$ and $r$ provide infinitely many
quadruples with the specific property, since $a$ has infinitely many choices,
by using Lemma 2.3 $(i)$.
Case II:  $2\mid m$ and $2\mid k$.
Here, we use $\alpha_{1}$ and $\alpha_{2}$ as defined in the previous case. We
consider $a=4(6a_{1}+1,6b_{1})$ with $Nm((6a_{1}+1,6b_{1}))=1$. Then (4.1)
gives

$r=((6Mm+6M+2m+(d/2)(6Nk+k)-12a_{1},6Mk+2k+3Nm+(m/2)-12b_{1}))$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{4}\Big{\\{}((6Mm+6M+2m+(d/2)(6Nk+k)-12a_{1})^{2}+d(6Mk+2k+3Nm+(m/2)-$
$\displaystyle
12a_{1})^{2}-4m-2,2(6Mm+6M+2m+(d/2)(6Nk+k)-12a_{1})(6Mk+2k+3Nm+$
$\displaystyle(m/2)-12b_{1})-4k)\times(6a_{1}+1,6b_{1})\Big{\\}}.$

Now for $r,b\in\mathbb{Z}[\sqrt{d}]$, we must have $m\equiv 2\pmod{4}$. Assume
that

$(\alpha,\beta)=(6Mm+6M+2m+(d/2)(6Nk+k),6Mk+2k+3Nm+m/2).$

Then $r=(\alpha-12a_{1},\beta-12b_{1})$. Now if $a+4b+4r=0$, then

$\displaystyle 4+\alpha^{2}+d\beta^{2}-4m-2+4\alpha=0,$ $\displaystyle
2\alpha\beta-4k+4\beta=0.$

By Lemma 2.4, we conclude that there exist only finitely many $\alpha$ and
$\beta$ which satisfy the above system of equations. We now rewrite,

$\displaystyle\alpha$ $\displaystyle=6M(m+1)+N(3dk)+2m+(d/2)k$
$\displaystyle\beta$ $\displaystyle=6Mk+3Nm+(m/2)+2k.$

These can be written as

$\begin{pmatrix}\alpha-2m-(d/2)k\\\
\beta-(m/2)-2k\end{pmatrix}=\begin{pmatrix}6(m+1)&3dk\\\
6k&3m\end{pmatrix}\begin{pmatrix}M\\\ N\end{pmatrix}.$

Since $m\equiv 2\pmod{4}$ and $k$ is even, so that the determinant of

$\begin{pmatrix}6(m+1)&3dk\\\ 6k&3m\end{pmatrix}$

is non-zero. As we have infinitely many choices of $M$ and $N$, so that there
exist infinitely many $\alpha$ and $\beta$ for which $a+4b+4r\neq 0$. Hence we
can take such $M$ and $N$ for which $a+4b+4r\neq 0$. Using these values of
$a,b$ and $r$, one can get infinitely many quadruples with the property
$D(n)$, since we have infinitely many choices of $a$, by using Lemma 2.3
$(i)$.

We now consider $m\equiv 0\pmod{4}$. In this case, we replace $n$ by
$n=(16m+2,8k)$, and we use the following factorization:

$3n=6(8m+1,4k)=(-12M-2,6N+1)(12M+2,6N+1)(8m+1,4k).$ (4.2)

Here, we choose

$\displaystyle\alpha_{1}$ $\displaystyle=(-12M-2,6N+1),$
$\displaystyle\alpha_{2}$
$\displaystyle=(96Mm+12M+16m+2+d(24Nk+4k),48Mk+8k+48Nm+8m+6N+1).$

Then

$a+2r=(48Mm+8m+d(12Nk+2k),24Mk+4k+24Nm+4m+6N+1).$ (4.3)

We take $a=(12a_{1}+4,6b_{1}+1)$ with $Nm(a)=6$, by Lemma 2.3 $(iv)$ then

$r=(24Mm+4m+d(6Nk+k)-6a_{1}-2,12Mk+2k+12Nm+2m+3N-3b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{6}\Big{\\{}((24Mm+4m+d(6Nk+k)-6a_{1}-2)^{2}+d(12Mk+2k+12Nm+2m+3N-$
$\displaystyle 3b_{1})^{2}-16m-2,2(24Mm+4m+d(6Nk+k)-6a_{1}-2)(12Mk+2k+12Nm+2m$
$\displaystyle+3N-3b_{1})-8k)\times(12a_{1}+4,-6b_{1}-1)\Big{\\}}.$

If $b,r\in\mathbb{Z}[\sqrt{d}]$, then

$(m,k)\equiv(0,1),(0,2),(1,0),(1,1),(2,0),(2,2)\pmod{(3,3)}.$

Let $a=(12a_{1}-4,6b_{1}+1)$ such that $Nm(a)=6$. Then (4.3) gives

$r=(24Mm+4m+d(6Nk+k)-6a_{1}+2,12Mk+2k+12Nm+2m+3N-3b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{6}\Big{\\{}((24Mm+4m+d(6Nk+k)-6a_{1}+2)^{2}+d(12Mk+2k+12Nm+2m+3N-3b_{1})^{2}$
$\displaystyle-16m-2,2(24Mm+4m+d(6Nk+k)-6a_{1}+2)(12Mk+2k+12Nm+2m+3N-$
$\displaystyle 3b_{1})-8k)\times(12a_{1}-4,-6b_{1}-1)\Big{\\}}.$

If $b,r\in\mathbb{Z}[\sqrt{d}]$, then we must have

$(m,k)\equiv(1,2)\pmod{(3,3)}.$

Again we utilize (4.2) by considering

$\begin{cases}\alpha_{1}=(12M+2,6N+1),\\\
\alpha_{2}=(-96Mm-12M-16m-2+d(24Nk+4k),-48Mk-8k+48mN+8m+6N+1).\end{cases}$

Then

$a+2r=(-48Mm-8m+d(12Nk+2k),-24Mk-4k+24mN+4m+6N+1).$ (4.4)

Here we choose $a=(12a_{1}+4,6b_{1}-1)$ such that $Nm(a)=6$. Then

$r=(-24Mm-4m+d(6Nk+k)-6a_{1}-2,-12Mk-2k+12mN+2m+3N+1-3b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{6}\Big{\\{}((-24Mm-4m+d(6Nk+k)-6a_{1}-2)^{2}+d(-12Mk-2k+12mN+2m+3N$
$\displaystyle+1-3b_{1})^{2}-16m-2,2(-24Mm-4m+d(6Nk+k)-6a_{1}-2)(-12Mk-2k+$
$\displaystyle 12mN+2m+3N+1-3b_{1})-8k)\times(12a_{1}+4,-6b_{1}+1)\Big{\\}}.$

Now $b,r$ to be in $\mathbb{Z}[\sqrt{d}]$ only when

$(m,k)\equiv(2,1)\pmod{(3,3)}.$

Case III:  $2\nmid m$ and $2\mid k$.
Here we consider the following factorization:

$3n=(-12M-2,6N+1)(12M+2,6N+1)(2m+1,2k).$

Let

$\displaystyle\alpha_{1}$ $\displaystyle=(-12M-2,6N+1),$
$\displaystyle\alpha_{2}$
$\displaystyle=(24Mm+12M+4m+2+d(12Nk+2k),24Mk+4k+12Nm+2m+6N+1).$

Then we have

$a+2r=(12Mm+2m+d(6Nk+k),12Mk+2k+6Nm+m+6N+1).$ (4.5)

We choose $a=2(6a_{1}+1,6b_{1})$ such that the norm of $(6a_{1}+1,6b_{1})$ is
one, by Lemma 2.3 $(i)$. Then

$r=(6Mm+m+(d/2)(6Nk+k)-6a_{1}-1,6Mk+k+3Nm+(m+1)/2+3N-6b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{2}\Big{\\{}((6Mm+m+(d/2)(6Nk+k)-6a_{1}-1)^{2}+d(6Mk+k+3Nm+(m+1)/2+$
$\displaystyle 3N-6b_{1})^{2}-4m-2,2(6Mm+m+(d/2)(6Nk+k)-6a_{1}-1)(6Mk+k+3Nm$
$\displaystyle+(m+1)/2+3N-6b_{1})-4k)\times(6a_{1}+1,-6b_{1})\Big{\\}}.$

Since $2\nmid m$ and $2\mid k$, so that $r,b\in\mathbb{Z}[\sqrt{d}]$.
Case IV:  Both $2\nmid km$.
We use $\alpha_{1}$ and $\alpha_{2}$ as defined in Case III. Here, we choose
$a=4(6a_{1}+1,6b_{1})$ with $Nm((6a_{1}+1,6b_{1}))=1$, by Lemma 2.3 $(i)$. Now
utilizing (4.5), we get

$r=(6Mm+m+(d/2)(6Nk+k)-12a_{1}-2,6Mk+k+3Nm+(m+1)/2+3N-12b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{4}\Big{\\{}((6Mm+m+(d/2)(6Nk+k)-12a_{1}-2)^{2}+d(6Mk+k+3Nm+(m+1)/2+$
$\displaystyle 3N-12b_{1})^{2}-4m-2,2(6Mm+m+(d/2)(6Nk+k)-12a_{1}-2)(6Mk+k+3Nm$
$\displaystyle+(m+1)/2+3N-12b_{1})-4k)\cdot(6a_{1}+1,-6b_{1})\Big{\\}}.$

If $r,b\in\mathbb{Z}[\sqrt{d}]$, then $m\equiv 3\pmod{4}$. The existence of
infinitely many quadruples can be seen by similar argument of $n=(4m+2,4k)$
with $m\equiv 2\pmod{4}$ and even $k$.

The next case is $m\equiv 1\pmod{4}$. In this case $n$ can be replaced by
$n=(16m+6,8k+4)$. Here we consider the following:

$3n=6(8m+3,4k+2)=(12M+4,-6N-1)(12M+4,6N+1)(8m+3,4k+2).$ (4.6)

Now we choose

$\displaystyle\alpha_{1}=$ $\displaystyle(12M+4,-6N-1),$
$\displaystyle\alpha_{2}=$
$\displaystyle(96Mm+36M+32m+12+d(24Nk+12N+4k+2),48Mk+24M+16k$
$\displaystyle+11+48Nm+18N+8m).$

Then we have

$\displaystyle a+2r=$
$\displaystyle(48Mm+24M+16m+8+d(12Nk+6N+2k+1),24Mk+12M+8k+5$
$\displaystyle+24Nm+6N+4m).$ (4.7)

Set $a=(12a_{1}+2,6b_{1}+1)$ with $Nm(a)=-6$, by Lemma 2.3 $(v)$. Then

$\displaystyle r=$
$\displaystyle(24Mm+12M+8m+4+(d/2)(12Nk+6N+2k+1)-6a_{1}-1,12Mk+6M+4k$
$\displaystyle+2+12Nm+3N+2m-3b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{-6}\Big{\\{}((24Mm+12M+8m+4+(d/2)(12Nk+6N+2k+1)-6a_{1}-1)^{2}+d(12Mk$
$\displaystyle+6M+4k+2+12Nm+3N+2m-3b_{1})^{2}-16m-6,2(24Mm+12M+8m+4+$
$\displaystyle(d/2)(12Nk+6N+2k+1)-6a_{1}-1)(12Mk+6M+4k+2+12Nm+3N+2m-$
$\displaystyle 3b_{1})-8k-4)\cdot(12M+2,-6N-1)\Big{\\}}.$

If $r,b\in\mathbb{Z}[\sqrt{d}]$, then

$(m,k)\equiv(0,0),(0,1),(1,1),(1,2),(2,0),(2,2)\pmod{(3,3)}.$

Let us choose $a=(12a_{1}-2,6b_{1}+1)$ with $Nm(a)=-6$, then (4) gives

$\displaystyle r=$
$\displaystyle(24Mm+12M+8m+4+(d/2)(12Nk+6N+2k+1)-6a_{1}+1,12Mk+6M+4k$
$\displaystyle+2+12Nm+3N+2m-3b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{-6}\Big{\\{}(((24Mm+12M+8m+4+(d/2)(12Nk+6N+2k+1)-6a_{1}+1)^{2}+d(12Mk$
$\displaystyle+6M+4k+2+12Nm+3N+2m-3b_{1})^{2}-16m-6,2((24Mm+12M+8m+4+$
$\displaystyle(d/2)(12Nk+6N+2k+1)-6a_{1}+1)(12Mk+6M+4k+2+12Nm+3N+2m-$
$\displaystyle 3b_{1})-8k-4)\cdot(12M-2,-6N-1)\Big{\\}}.$

Once again if $b,r\in\mathbb{Z}[\sqrt{d}]$, then
$(m,k)\equiv(0,2)\pmod{(3,3)}$.

Using (4.6), we again choose

$\displaystyle\alpha_{1}=$ $\displaystyle(12M+4,6N+1),$
$\displaystyle\alpha_{2}=$
$\displaystyle(96Mm+36M+32m+12+d(-24Nk-12N-4k-2),48Mk+24M+16k+8-$
$\displaystyle 48Nm-18N-8m-3).$

Then we get

$\displaystyle a+2r=$
$\displaystyle(48Mm+24M+16m+8+d(-12Nk-6N-2k-1),24Mk+12M+8k+4-$ $\displaystyle
24Nm-6N-4m-1).$

Here we set $a=(12a_{1}+2,6b_{1}-1)$ with $Nm(a)=-6$, by Lemma 2.3 $(v)$ then

$\displaystyle r=$
$\displaystyle(24Mm+12M+8m+4+(d/2)(-12Nk-6N-2k-1)-6a_{1}-1,12Mk+6M+$
$\displaystyle 4k+2-12Nm-3N-2m-3b_{1})$

and

$\displaystyle b=$
$\displaystyle\dfrac{1}{-6}\Big{\\{}((24Mm+12M+8m+4+(d/2)(-12Nk-6N-2k-1)-6a_{1}-1)^{2}+d(12Mk$
$\displaystyle+6M+4k+2-12Nm-3N-2m-3b_{1})^{2}-16m-6,2(24Mm+12M+8m+4+$
$\displaystyle(d/2)(-12Nk-6N-2k-1)-6a_{1}-1)(12Mk+6M+4k+2-12Nm-3N-2m-$
$\displaystyle 3b_{1})-8k-4)\cdot(12M+2,-6N+1)\Big{\\}}.$

For $b,r\in\mathbb{Z}[\sqrt{d}]$, we have $(m,k)\equiv(1,0)\pmod{(3,3)}$. This
completes the proof of Theorem 1.2.

##  5\. Concluding Remarks

Given a square-free integer, the existence of Diophantine quadruple in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$ for certain $n^{\prime}s$ has
been investigated in [6, 17]. We investigate the above problem for the
remaining values of $n$. However, our method does not work for a few values of
$n$, viz. $n=\\{4(12r+5,6s+3),4(12r+11,6s+3),(48r+38,24s+12),(48r+2,24s)\\}$
with $r,s\in\mathbb{Z}$.

We look for the existence of Diophantine quadruple in $\mathbb{Z}[\sqrt{d}]$
with the property $D(n)$ for the above exceptions. We first shorten these
exceptions with the help of [6, Theorem 1.1], and then we provide some
examples for the remaining cases. Let $d=2N$ with $N\in\mathbb{N}$ such that
equations (1.1) and (1.2) are solvable in integers. Assume that
$n=4(12m+5,6k+3)$ with $m=\alpha N+\beta$ and $k=\alpha_{1}N+\beta_{1}$, where
$\alpha,\beta,\alpha_{1},\beta_{1}\in\mathbb{Z}$. Then $n=(12\alpha
N+12\beta+5,6\alpha_{1}N+6\beta_{1}+3)$. Utilizing (iii) of Lemma 2.3, we get
$2,3\nmid N$ and thus we can choose $\beta,\beta_{1}$ such that $12\beta+5$
and $6\beta_{1}+3$ are of the form $N\gamma$ and $N\gamma_{1}$, respectively
with odd integer $\gamma$ and $\gamma_{1}$. Thus
$n=2N(24\alpha+2\gamma,12\alpha_{1}+2\gamma_{1})$. Since
$2\gamma,2\gamma_{1}\equiv 2\pmod{4}$, so that $24\alpha+2\gamma$ and
$12\alpha_{1}+2\gamma_{1}$ are of the form $4t_{1}+2$. Again $2N$ is square in
$\mathbb{Z}[\sqrt{d}]$ and thus [6, Theorem 1.1] and Lemma 2.1 together show
that there exist infinitely many Diophantine quadruples in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$. Analogously, one can draw a
similar conclusion for $n=4(12m+11,6k+3)$. We now consider
$n=(4(12m+9)+2,4(6k+3))$. As in the above,
$n=2(24N\alpha+24\beta+19,12N\alpha_{1}+12\beta_{1}+6)$. Since $2,3\nmid N$,
so we can choose $\beta,\beta_{1}$ such that $24\beta+19$ and $12\beta_{1}+6$
are of the form $N\gamma$ and $N\gamma_{1}$, respectively. Using (iii) of
Lemma 2.3, we get $N\equiv 1\pmod{4}$ and thus $\gamma\equiv 3\pmod{4}$ and
$\gamma_{1}\equiv 2\pmod{4}$. Finally we use [6, Theorem 1.1] and Lemma 2.1 to
conclude that there exist infinitely many Diophantine quadruples in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$. In this similar fashion, we
can establish a similar result for $n=(4(12m)+2,4(6k))$. We now provide some
examples supporting the existence of Diophantine quadruple in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$ for the remaining values of
$n$.

### Example 1

We consider $d=10$ and $n=4(12m+5,6k+3)$ with $m,k\in\mathbb{Z}$. Take $m=5M$
and $k=5K+2$, where $M,K\in\mathbb{Z}$. Then $n=4(5(12M+1),30k+15)$ which can
be written as $n=10(24M+2,12K+6)$. Thus $n$ is of the form
$4(4m^{\prime}+2,4k^{\prime}+2)$ with $m^{\prime},k^{\prime}\in\mathbb{Z}$.
Therefore using [6, Theorem 1.1] and Lemma 2.1, we conclude that there exist
infinitely many Diophantine quadruples in $\mathbb{Z}[\sqrt{10}]$ with the
property $D(n)$. Analogously, we can show the same for $n=4(12m+11,6k+3)$ by
putting $m=5M+2$ and $k=5K+2$.

### Example 2

Suppose $d=10$ and $n=(4(12m+9)+2,4(6k+3))=2(24m+19,12k+6)$. Take $m=5M+4$ and
$k=5K+2$. Then $n=10(24M+23,12k+6)$. Since $24M+23\equiv 3\pmod{4}$ and
$12K+6\equiv 2\pmod{4}$, so that by [6, Theorem 1.1] and Lemma 2.1, we can
conclude that there exist infinitely many Diophantine quadruples in
$\mathbb{Z}[\sqrt{10}]$ with the property $D(n)$. Analogously, we can show the
same for $n=(4(12m)+2,4(6k))$ by taking $m\equiv 1\pmod{5}$ and $k\equiv
0\pmod{5}$.

### Example 3

Assume that $d=10$ and $n=4(12m+5,6k+3)$ with $m\equiv 2,3\pmod{5}$. We
factorize

$\displaystyle 3n$ $\displaystyle=12(12m+5,6k+3)$
$\displaystyle=(-18,6)(3,1)(24m+10,12k+6)$
$\displaystyle=(-18,6)(120k+72m+90,~{}36k+24m+28).$

Now, we take $\alpha_{1}$ and $\alpha_{2}$ to the first and second factor of
the above equation, respectively. We utilize Lemma 2.2 to get

$a+2r=(60k+36m+36,18k+12m+17).$

We choose $a=(19,6)^{n}(0,1)$ with norm $-10$, where $n\in\mathbb{N}$. This
implies that there exist $\alpha,\beta\in\mathbb{Z}$ such that
$a=(20\alpha,10\beta-1)$. and thus
$r=(30k+18m-10\alpha+18,~{}9k+6m+9-5\beta)$. The equation $ab+n=r^{2}$ implies

$b=\dfrac{(r^{2}-n)(10\alpha,-10\beta-1)}{-10}.$

Since $m\equiv 2,\text{~{}~{}or~{}~{}}3\pmod{5}$, $b\in\mathbb{Z}[\sqrt{10}]$
and we have infinitely many $a$, therefore by using Lemmas 2.2 and 2.4, we get
infinitely many Diophantine quadruples in $\mathbb{Z}[\sqrt{10}]$ with the
property $D(n)$. Analogously, we can show the existence of infinitely many
Diophantine quadruples in $\mathbb{Z}[\sqrt{10}]$ with the property $D(n)$ for
$n=4(12m+11,6k+3)$ when $m\equiv 0,\text{~{}~{}or~{}~{}}4\pmod{5}$.

### Example 4

Suppose that $d=10$ and $(4(12m+9)+2,~{}4(6k+3))$ with $m\equiv
1,\text{~{}~{}or~{}~{}}2\pmod{5}$. We factorize

$\displaystyle 3n$ $\displaystyle=(4,1)(4,-1)(2(12m+9)+1,~{}2(6k+3))$
$\displaystyle=(4,1)(-120k+96m+16,~{}48k-24m+5).$

We choose $\alpha_{1}$ and $\alpha_{2}$ to the first and second factor of the
last equation, respectively. We use Lemma 2.2 to get
$a+2r=(-60k+48m+10,~{}24k-12m+3)$. Choose $a=(19,6)^{n}(10,3)$ with norm $10$,
where $n\in\mathbb{N}$. Thus there exist $\alpha,~{}\beta\in\mathbb{Z}$ such
that $a=(20\alpha+10,10\beta+3)$. So $r=(24m-30k-10\alpha,-6m+12k-5\beta)$.
Therefore Lemma 2.2 gives

$b=\dfrac{r^{2}-n}{a}.$

Since $m\equiv 1,\text{~{}~{}or~{}~{}}2\pmod{5}$, so that
$b\in\mathbb{Z}[\sqrt{10}]$. Hence there exist infinitely many Diophantine
quadruples in $\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$. Analogously, we
can deal quadruples with the property $D(n)$ for $n=(4(12m)+2,4(6k))$ when
$m\equiv 3,\text{~{}~{}or~{}~{}}4\pmod{5}$.

The problem of existence of infinitely many Diophantine quadruples in
$\mathbb{Z}[\sqrt{d}]$ with the property $D(n)$ for $n\in\mathbb{Z}[\sqrt{d}]$
is solved except for $n\in\mathcal{S}_{0}:=\\{S_{1},S_{2},S_{3},S_{4}\\}$,
where
$S_{1}=\\{4(12m+5,6k+3):(m,k)\equiv(0,0),(0,1),(0,3),(0,4)\pmod{(5,5)}\text{
and }m\equiv 1,4\pmod{5}\\}$,
$S_{2}=\\{4(12m+11,6k+3):(m,k)\equiv(2,0),(2,1),(2,3),(2,4)\pmod{(5,5)}\text{
and }m\equiv 1,3\pmod{5}\\}$,
$S_{3}=\\{(4(12m+9)+2,4(6k+3)):(m,k)\equiv(4,0),(4,1),(4,3),(4,4)\pmod{(5,5)}\text{
and }m\equiv 0,3\pmod{5}\\}$ and
$S_{4}=\\{48m+2,36k):(m,k)\equiv(1,1),(1,2),(1,3),(1,4)\pmod{(5,5)}\text{ and
}m\equiv 0,2\pmod{5}\\}$.

Finally we put the following question for the remaining values of $n$.

######  Question 5.1.

Are there infinitely many Diophantine quadruples in $\mathbb{Z}[\sqrt{10}]$
with the property $D(n)$ for $n\in S_{1}\cup S_{2}\cup S_{3}\cup S_{4}$?

## Acknowledgements

A. Hoque acknowledges SERB MATRICS Project (MTR/2021/000762), Govt. of India.

## References

  * [1] N. Adžaga, On the size of Diophantine $m$-tuples in imaginary quadratic number rings, Bull. Math. Sci. 9 (2019), no. 3, Article ID: 1950020, 10pp.
  * [2] F. S. Abu Muriefah and A. Al Rashed, Some Diophantine quadruples in the ring $\mathbb{Z}[\sqrt{-2}]$, Math. Commun. 9 (2004), 1–8.
  * [3] A. Baker and H. Davenport, The equations $3x^{2}-2=y^{2}$ and $8x^{2}-7=z^{2}$, Quart. J. Math. Oxford Ser. (2) 20 (1969), 129–137.
  * [4] M. Bliznac Trebješanin and A. Filipin, Nonexistence of $D(4)$-quintuples, J. Number Theory 194 (2019), 170–217.
  * [5] E. Brown, Sets in which $xy+k$ is always a square, Math. Comp. 45 (1985), 613–620.
  * [6] K. Chakraborty, S. Gupta, A. Hoque, On a conjecture of Franušić and Jadrijević: Counter-examples, Results Math. 78 (2023), no. 1, 14pp, article no. 18.
  * [7] K. Chakraborty, S. Gupta and A. Hoque, Diophantine triples with the property $D(n)$ for distinct $n$’s, Mediterr. J. Math. 20 (2023), no. 1, 13pp, article no. 31.
  * [8] A. Dujella, Generalization of a problem of Diophantus, Acta Arith. 65 (1993), 15–27.
  * [9] A. Dujella, The problem of Diophantus and Davenport for Gaussian integers, Glas. Mat. Ser. III 32 (1997), 1–10.
  * [10] A. Dujella, There are only finitely many Diophantine quintuples, J. Reine Angew. Math. 566 (2004), 183–214.
  * [11] A. Dujella, Number Theory, Školska knjiga, Zagreb, 2021.
  * [12] A. Dujella, Triples, quadruples and quintuples which are $D(n)$-sets for several $n$’s, in Class Groups of Number Fields and Related Topics, eds. K. Chakraborty, A. Hoque and P. P. Pandey (to appear).
  * [13] C. Elsholz, A. Filipin and Y. Fujita, On Diophantine quintuples and $D(-1)$-quadruples, Monats. Math. 175 (2014), 227–239.
  * [14] S. Gupta, $D(-1)$ tuples in imaginary quadratic fields, Acta Math. Hungar. 164 (2021), 556–569.
  * [15] B. He, A. Togbé and V. Ziegler, There is no Diophantine quintuple, Trans. Amer. Math. Soc. 371 (2019), 6665–6709.
  * [16] N. C. Bonciocat, M. Cipu and M. Mignotte, There is no Diophantine $D(-1)$-quadruple, J. London Math. Soc. 105 (2022), 63–99.
  * [17] Z. Franušić, Diophantine quadruples in the ring $\mathbb{Z}[\sqrt{2}]$, Math. Commun. 9 (2004), 141–148.
  * [18] Z. Franušić, Diophantine quadruples in $\mathbb{Z}[\sqrt{4k+3}]$, Ramanujan J. 17 (2008), 77–88.
  * [19] Z. Franušić, A Diophantine problem in $\mathbb{Z}[\sqrt{(1+d)/2}]$, Studia Sci. Math. Hungar. 46 (2009), 103–112.
  * [20] Z. Franušić, Diophantine quadruples in the ring of integers of the pure cubic field $\mathbb{Q}(\sqrt[3]{2})$, Miskolc Math. Notes 14 (2013), 893–903.
  * [21] Z. Franušić and B. Jadrijević, $D(n)$-quadruples in the ring of integers of $\mathbb{Q}(\sqrt{2},\sqrt{3})$, Math. Slovaca 69 (2019), 1263–1278.
  * [22] Z. Franušić and I. Soldo, The problem of Diophantus for integers of $\mathbb{Q}(\sqrt{-3})$, Rad Hrvat. Akad. Znan. Umjet. Mat. Znan. 18 (2014), 15–25.
  * [23] L. J. Matić, Non-existence of certain Diophantine quadruples in rings of integers of pure cubic fields, Proc. Japan Acad. Ser. A Math. Sci. 88 (2012), no. 10, 163–167.
  * [24] I. Soldo, On the existence of Diophantine quadruples in $\mathbb{Z}[\sqrt{-2}]$, Miskolc Math. Notes 14 (2013), 265–277.