Added Chi3Medium

CHANGELOG.md

@@ -7,6 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 - Added `hlim` and `vlim` kwargs to `Simulation.plot()` and `Simulation.plot_eps()` for setting horizontal and veritcal plot limits.
+- Added `Chi3Medium` for materials with a chi3 nonlinearity.
 
 ### Changed
 - `nyquist_step` also taking the frequency range of frequency-domain monitors into account.

tests/sims/simulation_2_4_0rc2.json

@@ -591,6 +591,33 @@
                 "subpixel": false
             }
         },
+        {
+            "geometry": {
+                "type": "Box",
+                "center": [
+                    -1.0,
+                    0.5,
+                    0.5
+                ],
+                "size": [
+                    1.0,
+                    1.0,
+                    1.0
+                ]
+            },
+            "name": null,
+            "type": "Structure",
+            "medium": {
+                "name": null,
+                "frequency_range": null,
+                "allow_gain": false,
+                "type": "Chi3Medium",
+                "permittivity": 1.0,
+                "conductivity": 0.0,
+                "numiters": 20,
+                "chi3": 0.1
+            }
+        },
         {
             "geometry": {
                 "type": "PolySlab",
@@ -1821,4 +1848,4 @@
     "normalize_index": 0,
     "courant": 0.8,
     "version": "2.4.0rc2"
-}
+}

tests/test_components/test_medium.py

@@ -508,3 +508,7 @@ def test_perturbation_medium():
             poles=[(1j, 3), (2j, 4)],
             poles_perturbation=[(None, pp_real)],
         )
+
+
+def test_nonlinear_medium():
+    med = td.Chi3Medium(chi3=1.5, numiters=20)

tests/utils.py

@@ -194,6 +194,13 @@
             ),
             medium=custom_sellmeier,
         ),
+        td.Structure(
+            geometry=td.Box(
+                size=(1, 1, 1),
+                center=(-1.0, 0.5, 0.5),
+            ),
+            medium=td.Chi3Medium(chi3=0.1, numiters=20),
+        ),
         td.Structure(
             geometry=td.PolySlab(
                 vertices=[(-1.5, -1.5), (-0.5, -1.5), (-0.5, -0.5)], slab_bounds=[-1, 1]

tidy3d/__init__.py

@@ -14,6 +14,7 @@
 from .components.medium import CustomMedium, CustomPoleResidue
 from .components.medium import CustomSellmeier, FullyAnisotropicMedium
 from .components.medium import CustomLorentz, CustomDrude, CustomDebye, CustomAnisotropicMedium
+from .components.medium import Chi3Medium
 from .components.transformation import RotationAroundAxis
 from .components.medium import PerturbationMedium, PerturbationPoleResidue
 from .components.parameter_perturbation import ParameterPerturbation
@@ -87,7 +88,7 @@
 from .material_library.parametric_materials import Graphene
 
 # for docs
-from .components.medium import AbstractMedium
+from .components.medium import AbstractMedium, AbstractNonlinearMedium
 from .components.geometry import Geometry
 from .components.source import Source, SourceTime
 from .components.monitor import Monitor

tidy3d/components/medium.py

@@ -34,6 +34,10 @@
 # extrapolation option in custom medium
 FILL_VALUE = "extrapolate"
 
+# max num iters for nonlinear medium
+# None for no max
+NONLINEAR_MAX_NUMITERS = None
+
 
 def ensure_freq_in_range(eps_model: Callable[[float], complex]) -> Callable[[float], complex]:
     """Decorate ``eps_model`` to log warning if frequency supplied is out of bounds."""
@@ -3186,6 +3190,73 @@ def perturbed_copy(
         return CustomPoleResidue.parse_obj(new_dict)
 
 
+class AbstractNonlinearMedium(AbstractMedium, ABC):
+    """Abstract nonlinear medium.
+
+    Note
+    ----
+    The nonlinear constitutive relation is solved iteratively; it may not converge
+    for strong nonlinearities. Increasing `numiters` can help with convergence.
+
+    """
+
+    numiters: pd.PositiveInt = pd.Field(
+        1,
+        title="Number of iterations",
+        description="Number of iterations for solving nonlinear constitutive relation.",
+    )
+
+    @pd.validator("numiters", always=True)
+    def _numiters_allowed(cls, val):
+        """Check that numiters is less than NONLINEAR_MAX_NUMITERS."""
+        if NONLINEAR_MAX_NUMITERS is None:
+            return val
+        if val > NONLINEAR_MAX_NUMITERS:
+            raise ValidationError(
+                "'AbstractNonlinearMedium.numiters' cannot be greater than "
+                f"{NONLINEAR_MAX_NUMITERS}."
+            )
+        return val
+
+
+class Chi3Medium(AbstractNonlinearMedium, Medium):
+    """Chi3 medium described by a chi3 nonlinear susceptibility.
+
+    Note
+    ----
+    The instantaneous nonlinear polarization is given by
+    .. math::
+
+        P_{NL} = \\epsilon_0 \\chi_3 |E|^2 E
+
+    Note
+    ----
+    The nonlinear constitutive relation is solved iteratively; it may not converge
+    for strong nonlinearities. Increasing `numiters` can help with convergence.
+
+    Note
+    ----
+    For complex fields (e.g. when using Bloch boundary conditions), the nonlinearity
+    is applied separately to the real and imaginary parts, so that the above equation
+    holds when both E and :math:`P_{NL}` are replaced by their real or imaginary parts.
+    The nonlinearity is only applied to the real-valued fields since they are the
+    physical fields.
+
+    Note
+    ----
+    Different field components do not interact nonlinearly. For example,
+    when calculating :math:`P_{NL}_x`, we approximate :math:`|E|^2 \\approx |E_x|^2`.
+    This approximation is valid when the E field is predominantly polarized along one
+    of the x, y, or z axes.
+
+    Example
+    -------
+    >>> medium = Chi3Medium(chi3=1)
+    """
+
+    chi3: float = pd.Field(..., title="Chi3", description="Chi3 nonlinear susceptibility.")
+
+
 # types of mediums that can be used in Simulation and Structures
 
 
@@ -3208,6 +3279,7 @@ def perturbed_copy(
     CustomAnisotropicMedium,
     PerturbationMedium,
     PerturbationPoleResidue,
+    Chi3Medium,
 ]
 
 

tidy3d/components/simulation.py

@@ -22,6 +22,7 @@
 from .medium import Medium, MediumType, AbstractMedium, PECMedium
 from .medium import AbstractCustomMedium, Medium2D, MediumType3D
 from .medium import AnisotropicMedium, FullyAnisotropicMedium, AbstractPerturbationMedium
+from .medium import AbstractNonlinearMedium
 from .boundary import BoundarySpec, BlochBoundary, PECBoundary, PMCBoundary, Periodic
 from .boundary import PML, StablePML, Absorber, AbsorberSpec
 from .structure import Structure
@@ -2885,6 +2886,12 @@ def make_eps_data(coords: Coords):
                     red_coords = Coords(**dict(zip("xyz", coords_reduced)))
                     eps_structure = get_eps(structure=structure, frequency=freq, coords=red_coords)
 
+                    if isinstance(structure.medium, AbstractNonlinearMedium):
+                        consolidated_logger.warning(
+                            "Evaluating permittivity of a nonlinear "
+                            "medium ignores the nonlinearity."
+                        )
+
                     if isinstance(structure.geometry, TriangleMesh):
                         consolidated_logger.warning(
                             "Client-side permittivity of a 'TriangleMesh' may be "