xarray.align can return read-only arrays

What happened?

The attached code snippet raises "ValueError: output array is read-only" because the d1 array returned from xarray.align is read-only.

What did you expect to happen?

d1 should not be a read-only. d2 should have been added to it without issue.

Minimal Complete Verifiable Example

import numpy as np
import xarray as xr

d1 = xr.DataArray(np.array([]), coords={'dim': []})
d2 = xr.DataArray(np.array([1]), coords={'dim': ['x']})
d1 , d2 = xr.align(d1, d2, join='outer', fill_value=0, copy=True)
d1 += d2

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

Traceback (most recent call last):
  File "/home/atest.py", line 7, in <module>
    d1 += d2
  File "/some/where/lib/python3.10/site-packages/xarray/core/_typed_ops.py", line 328, in __iadd__
    return self._inplace_binary_op(other, operator.iadd)
  File "/some/where/lib/python3.10/site-packages/xarray/core/dataarray.py", line 4659, in _inplace_binary_op
    f(self.variable, other_variable)
  File "/some/where/lib/python3.10/site-packages/xarray/core/_typed_ops.py", line 680, in __iadd__
    return self._inplace_binary_op(other, operator.iadd)
  File "/some/where/lib/python3.10/site-packages/xarray/core/variable.py", line 2710, in _inplace_binary_op
    self.values = f(self_data, other_data)
ValueError: output array is read-only

Anything else we need to know?

I suspect this is related to https://github.com/pydata/xarray/issues/3813, i.e. that align is calling expand_dims behind the scenes. But even if the expand_dims aspect is considered a feature rather than a bug, the DataArray.copy() required to fix it should be done by align since the documentation says it should return a copy, not a view.

Environment

Reproduce the MVCE with the reported xarray 2023.9.0 environment, then trace the aligned result through the reported xarray/core/dataarray.py, xarray/core/variable.py, and _typed_ops.py paths. Add a regression test showing that align(..., copy=True) returns writable data and that d1 += d2 succeeds.