Faraday Rotation and Depolarization
====================================

The Faraday rotation effect
---------------------------

When **linearly polarized** radio emission (Stokes Q and U) travels through **magneto-ionic gas** — ionized plasma threaded by a magnetic field — the electric vector rotates. The rotation angle depends on the square of the wavelength, :math:`\lambda^2`, and on the **rotation measure** (RM), which is the line-of-sight integral of the electron density times the parallel component of the magnetic field. So the observed polarization position angle changes with frequency in a predictable way. This effect is used to probe magnetic fields and thermal plasma in the interstellar medium, in galaxies, and in galaxy clusters.

What we model in simulation is exactly this: intrinsic (Q, U) from a polarized source, transformed by the magneto-ionic medium into observed (Q', U') as a function of :math:`\lambda^2`. Depending on how the emission and the medium are distributed in **Faraday depth** (i.e. in RM along the line of sight), we get either pure rotation or rotation plus **depolarization**.

.. figure:: diagrams/png/fig09_models_sky_faraday-1.png
   :alt: Sky source and Faraday model hierarchy
   :width: 100%

   **Sky/Faraday model hierarchy** (sources, Faraday components, and composite relationships).

Faraday thin vs Faraday thick
-----------------------------

- **Faraday thin:** The polarized emission effectively passes through a **single Faraday depth**. That can mean: (1) a discrete **Faraday screen** — a slab of magneto-ionic gas in front of the source that rotates the polarization by one effective RM; or (2) an **emitting region** that is so narrow in Faraday depth (e.g. a compact source behind a uniform screen) that all rays see the same rotation. In either case, the position angle rotates with :math:`\lambda^2` (you measure one RM), and the **polarized intensity does not drop** with wavelength; there is no depolarization.

- **Faraday thick:** The emission is associated with a **range of Faraday depths**. That can mean: (1) an **extended magneto-ionic medium** (e.g. a turbulent or stratified screen) where different lines of sight have different RMs; or (2) **internal Faraday dispersion** — the source itself is extended in Faraday depth (emission and rotation mixed). Then, at a given :math:`\lambda^2`, contributions from different Faraday depths add with different phases, so the net polarization **rotates and depolarizes** with wavelength. The observed signal is a blend over a Faraday depth distribution; the wider the distribution, the stronger the wavelength-dependent depolarization.

In Pyralysis, **ThinFaradayComponent** implements the thin case (single RM, rotation only), and **ThickFaradayComponent** implements the thick case (parametrized Faraday depth distribution). You can chain components with **FaradayComposite** (e.g. thin screen then thick screen). All components operate on Dask arrays and are applied during visibility simulation when you attach them to a polarized source via :meth:`add_faraday_component <pyralysis.models.sky.source.SkySource.add_faraday_component>`.

Thin component
--------------

**ThinFaradayComponent** implements the thin case: rotation by a single rotation measure :math:`\mathrm{RM}` (rad/m²):

.. math::
   Q' = Q \cos(2\,\mathrm{RM}\,\lambda^2) - U \sin(2\,\mathrm{RM}\,\lambda^2)
   \qquad
   U' = Q \sin(2\,\mathrm{RM}\,\lambda^2) + U \cos(2\,\mathrm{RM}\,\lambda^2)

Polarized intensity is preserved.

**Parameters:**

- **rm** (float): Rotation measure in rad/m².

**Example:**

.. code-block:: python

   from pyralysis.models.faraday import ThinFaradayComponent
   comp = ThinFaradayComponent(rm=25.0)
   # Use with a polarized source: source.add_faraday_component(comp)

Thick component
--------------

**ThickFaradayComponent** implements the thick case. The complex response :math:`F(\lambda^2)` multiplies the intrinsic polarization :math:`P = Q + iU` so that :math:`P' = P \cdot F`.

**Parameters:**

- **rm_mean** (float): Mean Faraday depth in rad/m² (rotation and center of distribution).
- **rm_width** (float): Width of the distribution in rad/m² (must be non-negative).
- **kind** (``"sinc"`` or ``"gaussian"``, optional): Shape of the Faraday depth distribution. Default is ``"sinc"``.

**Sinc / top-hat** (``kind="sinc"``, default): Top-hat Faraday depth distribution over
:math:`[\mu - \delta, \mu + \delta]` rad/m² with :math:`\delta` = **rm_width**. Response:

.. math::
   F = \exp(2i\mu\lambda^2) \,\mathrm{sinc}(2\delta\lambda^2)

**Gaussian** (``kind="gaussian"``): Faraday depth distribution
:math:`\phi(\mathrm{RM}) \propto \exp(-(\mathrm{RM} - \mu)^2 / (2\sigma^2))` with :math:`\sigma` = **rm_width**. Response:

.. math::
   F = \exp(2i\mu\lambda^2) \exp(-2\sigma^2\lambda^4)

**Example:**

.. code-block:: python

   from pyralysis.models.faraday import ThickFaradayComponent
   # Default: top-hat (sinc), half-width 10 rad/m²
   thick_sinc = ThickFaradayComponent(rm_mean=50.0, rm_width=10.0)
   # Gaussian distribution, mean RM 50 rad/m², sigma 10 rad/m²
   thick_gauss = ThickFaradayComponent(rm_mean=50.0, rm_width=10.0, kind="gaussian")

Composing components
-------------------

Use the ``+`` operator to apply multiple components in sequence. The output (Q', U') of one component becomes the input (Q, U) of the next. You can also use ``-`` to compose with the *inverse* of a component (e.g. ``thin1 + (-thin2)`` or ``thin1 - thin2``).

.. code-block:: python

   from pyralysis.models.faraday import ThinFaradayComponent, ThickFaradayComponent
   # Thin screen (RM 25) then thick (default sinc; mean 50, width 10)
   faraday = ThinFaradayComponent(25.0) + ThickFaradayComponent(50.0, 10.0)
   source.add_faraday_component(faraday)

Use in simulation
-----------------

1. Define a **polarized** sky source. There are two ways:

   - **Stokes vector:** :class:`PointSource <pyralysis.models.sky.point_source.PointSource>` with ``reference_intensity=[I, Q, U, V]`` (Stokes parameters in Jy).
   - **Polarization fraction and angle:** :meth:`PointSource.from_polarization <pyralysis.models.sky.source.Source.from_polarization>` with ``total_intensity``, ``pol_fraction``, ``pol_angle`` (and optional ``circular_pol_fraction``); pass ``sky_position``, ``reference_frequency``, etc. in ``**kwargs``.

2. Create one or more Faraday components and compose them with ``+`` if needed.
3. Attach the (composite) component to the source with :meth:`add_faraday_component <pyralysis.models.sky.source.SkySource.add_faraday_component>`.
4. Run the simulator as usual; the component(s) are applied when generating visibilities.

**Example (Stokes vector):**

.. code-block:: python

   from pyralysis.models.sky import PointSource
   from pyralysis.models.faraday import ThinFaradayComponent, ThickFaradayComponent
   from pyralysis.simulation import Simulator

   source = PointSource(
       reference_intensity=[1.0, 0.3, 0.1, 0.0],  # [I, Q, U, V] in Jy
       sky_position="12:00:00 45:00:00",
       reference_frequency=1e9,
   )
   faraday = ThinFaradayComponent(25.0) + ThickFaradayComponent(50.0, 10.0)
   source.add_faraday_component(faraday)

   sim = Simulator(interferometer=interferometer, sources=source)
   dataset = sim.simulate(create_dataset=True)

**Example (polarization fraction and angle):**

.. code-block:: python

   from astropy import units as u
   from pyralysis.models.sky import PointSource
   from pyralysis.models.faraday import ThinFaradayComponent, ThickFaradayComponent

   source = PointSource.from_polarization(
       1.0 * u.Jy,
       pol_fraction=0.1,
       pol_angle=45 * u.deg,
       sky_position="12:00:00 45:00:00",
       reference_frequency=1e9,
   )
   faraday = ThinFaradayComponent(25.0) + ThickFaradayComponent(50.0, 10.0)
   source.add_faraday_component(faraday)

See also
--------
- :doc:`simulation` — Simulation framework and sky models
- :ref:`defining-sky-models` — Defining sky models (simulation.rst)
- API: :py:mod:`pyralysis.models.faraday`
